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A FINITE-CIRCUIT-ELEMENT  CODE 


FOR  MODELING  THE  COMPRESSION  OF  A GYRATING 
CHARGED-PARTICLE  BEAM 

I.  INTRODUCTION 


Over  the  .ast  several  years  a great  deal  of  interest  has  arisen  in  connection  with  the  topic 
of  gyrating  intense  ion  beams.11  11  A ring  or  cylindrical  current  layer  is  produced  by  the  motion 
of  the  ions  in  the  superposed  background  (quasiuniform)  magnetic  field  and  the  poloidal  self- 
field, with  ring  major  radius  R equal  to  the  ion  gyroradius.  If  the  net  current  in  such  a 


configuration  is  strong  enough,  the  direction  of  the  field  lines  within  the  ring  can  be  opposite 
that  of  the  background  field  (Fig.  1).  When  the  poloidal  field  on  axis,  Bp  = n„l /2R, 
exceeds  the  background  field  fl„,  the  field  in  the  interior  region  is  completely  reversed. 

Recently  it  has  been  proposed  to  increase  the  intensity  of  the  neutral  beams  used  to  heat 
the  plasma  in  2xllB  and  similar  mirror  devices  in  order  to  produce  field  reversal.14'  As  pointed 
out  by  Baldwin  and  Rensink,'5'  electric  fields  induced  by  the  buildup  of  current  tend  to  partially 
cancel  the  ion  current.  It  is  thus  unclear  that  an  initially  unreversed  configuration  can  become 
reversed,  no  matter  how  much  ion  current  is  added.  Even  if  the  configuration  is  compressed 
radially  (by  the  action,  e g.,  of  external  coils,  an  imploding  liner,  or  axial  translation  in  a tank 
with  converging  metal  >valls),  field  reversal  is  problematic.  The  flux  linking  the  ion  ring  tends 
to  be  conserved,  and  collisional  diffusion  only  flattens  the  profiles. 
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The  present  paper  describes  a code  developed  for  treating  the  dynamics  of  a gyrating  ion 
ring  interacting  with  a background  plasma  and  a (possibly  imploding)  metal  liner.  The  code  is 
called  1PICAC  (for  Ion  Beam-Plasma  Interaction  with  Cylindrical  Adiabatic  Compression).  It  is 
two-dimensional  (in  r,  z)  and  assumes  axisymmetry,  but  does  not  employ  finite-differences  on  a 
2D  grid  to  solve  the  dynamical  problem.  Instead,  each  portion  of  the  system  which  carries 
current  is  regarded  as  part  of  a circular  current  loop.  The  beam  is  one  such  loop,  the  liner  or 
wall  may  be  approximated  by  several  loops  side  by  side.  These  current  loops  are  coupled  by 
their  mutual  inductances,  and  the  dynamical  behavior  is  determined  through  solution  of  the  cir- 
cuit equations.  Thus  the  system  is  described  by  ordinary  differential  equations,  rather  then  the 
partial  differential  equations  of  the  usual  magnetohydrodynamic  treatment. 

The  principal  difficulty  in  this  approach  lies  in  determining  the  inductances.  These  change 
as  the  geometry  of  the  beam  and  liner  changes,  and  have  to  be  recalculated  at  every  timestep. 
Unless  some  approximation  is  invoked  to  simplify  them,  no  computational  advantage  results 
from  the  circuit  theory  technique.  Fortunately,  such  an  approximation  is  available  in  many 
charged-particle  ring  configurations  of  interest,  namely  that  of  large  aspect  ratio.  That  is,  the 
major  radius  R , of  the  j th  current  loop  is  taken  to  be  large  compared  with  its  minor  dimension 
and  the  separation  in  the  r—z  plane  between  it  and  any  other  loop.  It  is  not  necessary  but  is 
often  convenient  to  assume  that  resistance  and  current  are  distributed  uniformly  throughout  the 
r-z  cross  section  of  the  loop.  The  latter  may  be  of  arbitrary  shape,  but  is  usually  taken  to  be 
circular  or  rectangular. 

In  this  conception,  collisions  between  the  ion  beam  and  background  plasma  enter  as  a 
resistance  (and  possibly  an  Ohkawa  current161).  Plasma  energy  losses  by  radiation  and  convec- 
tion also  affect  the  beam  dynamics  through  the  inductances  and  the  resistance.  Consistent  with 
this  approach,  the  inertia  of  the  various  particle  species  is  ignored  (except  in  the  centrifugal 
force),  so  that  the  beam  and  plasma  remain  in  force  balance  with  the  wall  currents. 
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The  code  described  here  was  originally11  71  developed  for  an  ion-beam-plasma  interaction 
problem  related  to  but  distinct  from  that  of  producing  field  reversal.  We  started  with  a ring  of 
deuterium  (D)  ions  assuming  an  already  existing  field-reversed  geometry.  The  ring  was 
compressed  by  implosion  of  the  liner,  and  the  thermonuclear  energy  production  arising  from 
collisions  between  the  beam  ions  and  T or  He3  target  ions  in  the  background  plasma  was  stu- 
died. An  attempt  was  made  to  balance  the  components  of  the  system  so  that  the  collisiona! 
slowing-down  of  the  beam  ions  just  canceled  their  tendency  to  speed  up  because  of  angular 
momentum  conservation.  "Clamping"  the  beam  in  this  way  at  the  energy  for  which  the  beam- 
target  reaction  rate  peaks  ( — -150  keV  for  D-T  reactions)  maximizes  Q , the  ratio  of  the  yield  to 
the  sum  of  liner  and  plasma  energy.  It  was  found  however  that  even  with  optimized  parame- 
ters, Q was  limited  to  10%  or  less.  The  reason  was  that  the  energy  given  up  by  beam  ions  in 
collisions,  most  of  which  went  into  electron  heating,  caused  expansion  of  the  toroidal  beam- 
plasma  system  and  reduced  all  of  the  number  densities,  and  accordingly  reduced  the  beam- 
target  reaction  rate.  Presumably  Q would  increase  if  a method  were  found  to  cool  the  electrons 
and  recycle  their  thermal  energy. 

Some  results  from  this  earlier  work  will  be  displayed  for  purposes  of  illustration,  but  the 
method  is  much  more  general  in  applicability.  Instead  of  assuming  a preexisting  state  of  field 
reversal,  one  can  employ  the  code  to  study  its  origin  and  development  in  time.  This  problem 
will  not  however,  be  addressed  in  the  present  paper,  which  is  devoted  to  describing  the  code 
and  some  of  the  techniques  employed  in  its  implementation.  The  plan  of  the  paper  is  as  fol- 
lows. In  Section  II  we  derive  the  equations  of  the  circuit  theory  model  of  the  beam-plasma- 
liner  dynamical  system.  In  Section  111  we  discuss  isentropic  (lossless)  compression  of  an  ion 
ring  and  the  role  of  the  induced  electron  current  in  the  resultant  scaling.  Collisions  are 
described  in  Section  IV.  The  implementation  of  conduction,  particle  transport  processes  and 
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other  phenomena  is  discussed  in  Section  V,  and  an  example  is  described  in  Section  VI.  Our 
results  are  summarized  in  Section  VII.  A listing  of  the  code  is  given  in  an  Appendix. 


II.  LINER  MOTION  AND  EQUIVALENT  CIRCUIT  EQUATIONS 

It  is  natural  to  represent  the  ion  beam  (and  the  currents  carried  by  the  electrons  and  tar- 
get ion)  as  a current  loop.  It  is  equally  convenient,  though  perhaps  less  natural,  to  represent 
the  axial  current  profile  on  the  liner  (and  possibly  on  the  driver  coil)  as  a superposition  of  coax- 
ial current  loops.  Each  such  loop  constitutes  an  electrical  circuit  individually  coupled  to  each  of 
the  others,  and  contains  a self-inductance  and  a resistance  (arising  from  charged-particle 
encounters  in  the  case  of  the  ring).  The  circuit  elements  vary  in  time  as  the  geometry  charges. 

Thus  it  is  possible  to  calculate  the  implosion  dynamics  to  any  desired  degree  of  realism 
entirely  by  means  of  the  equivalent  circuit  equations.  This  representation  is,  in  fact,  a type  of 

"finite-element"  simulation.  The  minimum  number  of  such  circuits  required  to  describe 
electromagnetic  implosions  of  the  liner  is  one  each  for  the  driver,  liner  and  ring.  In  this  limit 

the  equivalent  circuit  is  that  shown  in  Fig.  2. 

The  circuit  equations  take  the  form 


where  j runs  over  all  current-carying  loops  in  the  system.  For  the  circuit  of  Fig.  2,  j 
(signifying  driver,  liner,  and  ring,  respectively).  The  flux  threading  the  j th  element  is 


(1) 

d.l,r 


(2) 

k 

where  is  the  inductance  coupling  circuits  j and  k , and  is  the  resistance  of  the  j th  circuit. 
Equation  (l)  describes  the  evolution  of  <t>,.  Given  a knowledge  of  the  >1),  and  the  induction 
coefficients  )!?,*,  Eq.  (2)  then  can  be  solved  for  the  It  by  matrix  inversion. 
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If  the  driver  is  static  and  energized  only  during  the  outermost  portion  of  the  cycle,  we  can 
make  an  additional  simplification  by  restricting  our  attention  to  times  when  the  liner  and  ring 
are  far  removed  from  the  driver  coil.  Then  j,  k take  on  only  the  values  l,r , and  there  are  just 
two  each  of  equations  (1)  and  (2).  The  numerical  results  described  and  plotted  below  were 
obtained  using  this  two-loop  circuit.  It  should  be  clear,  however,  that  most  of  the  discussion 
which  follows  is  independent  of  the  number  of  loops  employed.  We  have  experienced  no 
difficulty  in  implementing  versions  of  the  code  where  as  many  as  ten  loops  are  employed  to 
simulate  the  current  profile  in  the  liner.  It  appears  that  it  would  be  easy  to  generalize  the 
method  to  multiple  ion  rings  or  single  rings  with  multiple  constituent  current  filaments. 

The  coefficients  %Jk  are  very  easily  calculated.  Since  the  ring  deforms  freely,  it  tends  to 
evolve  so  as  to  maximize  its  self-inductance,  that  is,  toward  a circular  cross-section.  Moreover, 
one  wants  to  consider  configurations  where  ring  and  liner  are  close  together,  to  minimize  the 
volume  filled  with  magnetic  energy.  Thus  all  distances  separating  current-carrying  filaments  are 
small  compared  with  the  major  radii  R,  R,  (Fig.  1).  In  this  limit  the  self  and  mutual  induc- 
tances can  be  calculated  in  the  large-aspect-ratio  approximation  as 

£•  5n  = n0R  [In  (SR)  - 2 - InD],  (3) 

where  the  average  is  over  the  current-carrying  part  of  the  cross  section,  and  the  minor  diameter 

D satisfies  D « R. 

Using  (3)  we  find  that  the  self-inductance  of  the  ring  is  given  by 

= M„*  [In(8/?/r)  - 2 + 8]  (4) 

where  r is  the  ring  minor  radius,  and  8 depends  on  the  details  of  the  assumed  current  profile. 

For  example,  if  all  the  current  is  carried  in  a skin  located  at  the  minor  radius,  8 — 0;  if  the 

current  is  uniformly  distributed,  8 — 0.25;  and  if  the  ring  looks  like  a Bennett  pinch  in  cross- 

section,  8 — 0.5.  Similarly,  the  self-inductance  of  a liner  segment  is  approximately  (assuming 

the  current  is  carried  on  the  inner  surface) 
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*1=  *#  =*M.«/[ln  (&R/0  - 1/2]  (5) 

where  / is  the  length  of  the  segment,  assumed  much  larger  then  the  thickness,  and  R , is  the 


inside  radius;  and 


% I'  “ In  1 - nyj-  -1  (6) 

[(jr-/*,)a  + (1/2)*] ' 

- [(/?,  - R)l(ll 2)]  tan -‘[(//2)/(J?,  - *)]J 

More  important  than  the  exact  forms  of  (5)  and  (6)  (which  depend  on  the  cross  sections  as- 
sumed to  describe  the  liner)  is  the  fundamental  geometrical  requirement  M}r  ^ LrLh  with 
equality  holding  only  if  R = R Since  Eqs.  (4-6)  are  approximate,  this  inequality  must  be  en- 
forced by  means  of  an  explicit  interpolation;  otherwise,  the  ring  can  pass  right  through  the 
liner.  The  interpolation  formula  actually  used  is 

* /,  = tofr  + [(£,£,)>'*  - %,r]  [l  + [(*,  - *)/,]'].  (7) 

whet. Af'  \s  the  corrected  value  of  the  mutual  inductance.  The  dynamical  results  are  not  very 
sensitive  to  the  choice  of  p,  which  was  taken  to  be  10  in  the  numerical  calculation. 

As  is  well  known  from  electromagnetic  theory,  the  force  tending  to  change  any  coordinate 
9 on  which  an  inductive  coefficient  %Jk  depends  is  given  by 

_ , , 9 HJk 

Fjk  Ijh  qq  ■ (8) 

Employing  (8)  consistently  with  the  definitions  used  for  Hi,*  guarantees  conservation  of  total 
energy,  the  magnetic  portion  of  which  is 

^4^  ^ A1  jh  = ~ 2 'Jj^j  (9) 

J*  1 j 

Thus  in  carrying  out  numerical  calculations,  we  determine  the  total  force  of  the  ring  acting  on 
the  liner  according  to 
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F,  - - A I/, 


11- 
dR,  ' 


(10) 


where  the  summation  runs  -er  the  ring  and  all  segments  of  the  liner;  while  the  same  expres- 
sion with  opposite  sign  yields  the  force  with  which  the  liner  tends  to  hold  the  ring  in  place.  The 
liner  equation  of  motion  is  thus 

M[Rt  = F,  (]1) 

Similarly,  the  electromagnetic  force  acting  to  constrict  the  ring  is  given  by  Eq.  (8)  with 
9 - r: 


(12) 


Most  of  the  force  F,  comes  from  the  term  containing  %rr  = £r.  Because  of  the  use  of  the 
interpolation  formula,  Eq.  (7),  however,  there  is  a small  contribution  from  the  liner-ring 
mutual  inductances. 


III.  ISENTROPIC  COMPRESSION 

It  is  possible  to  develop  scaling  laws  in  terms  of  which  the  liner  motion  and  beam  and 
plasma  evolution  are  described  by  analytic  expressions,  provided  we  assume  the  absence  of 
both  fusion  reactions  and  loss  mechanisms.  This  model  is  not  a useful  starting  point  about 
which  to  perturb  to  describe  a realistic  reactor  design,  because  the  latter  is  quite  sensitive  to 
beam  slowing  and  the  heating  resulting  from  production  of  charged  fusion  reaction  products.  It 
is,  however,  valuable  in  describing  the  dynamics  in  the  absence  of  a target  plasma,  as  well  as 
guiding  us  in  developing  an  intuition  about  the  interdependence  of  various  parts  of  the  system. 

If  the  liner  is  represented  by  J,  distinct  current-carrying  segments,  there  are  J , + 1 fluxes 
and  J,  + 11  physical  variables.  In  our  numerical  calculations  we  usually  took  7,-1.  For  this 
case  the  12  physical  quantities  used  to  describe  a dynamical  state  of  the  system  are  the  fluxes 
d),  and  <Pr  linking  the  liner  and  ring,  respectively;  R and  R,;  the  ring  minor  ^dius  r,  the  total 
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numbers  of  beam  and  target  ions,  NB  and  .V7,  respectively;  the  beam,  target  and  electron  tem- 
peratures, Th , T,  and  Tc,  respectively;  and  the  mean  azimuthal  ion  drift  velocities  vfi  and  vr. 
To  proceed,  we  write  down  all  the  conservation  laws  that  are  available.  The  conserved  quanti- 
ties are  the  magnetic  flux  threading  the  j th  liner  segment 

<t>,  = I.  '’'A//  + = *>?.  03) 

/ 

and  that  threading  the  ring, 

- £,/,  + £ XriI,  = <D,°;  (14) 

i 

the  specific  angular  momentum  of  beam  ions. 


R\  B = R°v°,  (15) 

and  of  target  ions, 

R v r = R “v°;  (16) 

the  total  ion  numbers  for  each  species 

Hb  - 07) 

NT-W\  08) 

and  the  beam,  target  and  electron  entropy  functions: 


Here 


rByy = r^(  v°)y-\ 
ttv y~l  = Tf  (v°)y-\ 

T,v = 7-°(  v°)y~'. 

V = 2 n2Rr2  is  the  volume  of  the  beam/plasma  ring.  Superscripts  (°) 


09) 

(20) 

(21) 

indicate  an  initial  or 


a reference  state  of  the  system  (e.g.,  the  state  of  maximum  compression).  To  these  equations 
must  be  added  the  condition  of  force  balance  on  the  ring  in  the  direction  of  major  and  minor 
radius, 


0 


8JPr 

6 R 


.2 


, v a 1 dV  , HBmBv2B  NTmT\ } 

" 1 L h ~ZR  + ' 1 ” + P - + =~  + 


6R 


(22) 


NRL  MEMORANDUM  REPORT  3827 


and 


o = — / 2 + / T / 

° 2 dr  + A f 

respectively.  Here  p — fclAgTa  + NtTt  + A,  7,,]  K_1  is  the  internal  pressure  in  the  ring  ( k is 
the  Boltzmann  constant),  and  the  electron  number  is  obtained  from  the  condition  of  charge 


a K/  ^ a v 

-T—  + P -T-. 
dr  dr 


(23) 


neutrality, 


Ne  — Nb^b  + Nt^T’  (24) 

where  Za  is  the  charge  state  of  ion  species  a.  The  last  two  terms  in  eq.  (22)  are  the  centrifugal 

force  terms  derived  from  the  circulation  of  the  respective  species;  that  corresponding  to  the  tar- 
get ions  is  usually  negligible. 

Equations  (22)  and  (23)  have  been  derived  assuming  that  the  ring  inertia  is  negligible, 
i.e.,  that  the  ring  repositions  itself  instantaneously  in  response  to  any  change  in  the  position  of 
the  liner.  In  addition,  the  electron  mass  has  been  set  to  zero  systematically,  as  negligible  in 
comparison  with  those  of  the  ions.  The  ring  current  /,  satisfies 

A - lB+  h + (25) 

where 


and 


Ib 

h 


A’BeZflvfi 

2-irR 

NTeZf/T 

2nR 


(26) 

(27) 


Nee\t 
2 nR  ' 


(28) 


Equations  (1),  (11)  and  (13-23)  contribute  a set  of  12  + J , fundamental  algebraic  equations  in 
terms  of  the  12  + 7/  physical  quantities  defining  the  state.  [All  the  others  are  expressible  in 
terms  of  these  through  Eqs.  (2),  (4-6),  and  (24-28).]  Thus,  specifying  the  state  variables 
determines  the  evolution  of  the  system  completely.  We  rewrite  the  liner  force  equation  as 
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;£<*/*/)  = R,2R?(Rr2  - Ri~2)  + I y/,2  1^-  + 1,1,  ^ /(2p/./e,)J/ln(*,7/?,2) 

which  parametrizes  the  dynamical  history  in  terms  of  t.  Equation  (29)  is  derived  by  assuming 
conservation  of  the  liner  mass  M,  = lirpL  ( /?, 2 — R?)\  p is  the  (uniform)  liner  density,  L is 
the  overall  length,  and  /?,  is  the  outer  liner  radius. 

Let  us  assume  now  that  the  electron  current  tending  to  neutralize  1B  is  zero.  Then  by 
conservation  of  angular  momentum. 


1,  — 2tiR  ^dWd  + ^t^t)  “ yyyy  (NdR\ p + NtRvt)  ~~  R 2- 


The  minor  radius  force  balance  condition  (23)  reduces  to 


P~T~V 


Equations  (19-21),  weighted  by  the  respective  total  numbers  Nj,  sum  to  the  adiabatic  law 


pVy  = const. 


Taking  y = 5/3  and  combining  (31)  and  (32)  yields 


Hence  the  number  densities  for  species  a(a  = B,  T,  e),  na  = NJV , satisfy 

na  ~ — R~9/2,  (34) 

and  the  poloidal  field  near  the  ring  Bv  = p„Ir/2nr  satisfies 

Bj  - p ~ R-'5'2  (35) 

We  thus  have  a situation  in  which  almost  three-dimensional  compression  of  the  ring  occurs  as 

R = Ri  is  reduced.  The  poloidal  field  (35)  rises  almost  as  the  inverse  fourth  power  of  R and 

the  temperatures  scale  like  T ~ /?'3. 

At  the  other  extreme,  the  motion  of  the  liner  may  be  such  as  to  induce  electron  currents 
le  totally  neutralizing  the  change  in  ion  current, 

/,  ~ const  (36) 
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Going  through  the  same  steps  as  above,  we  find 

r ~ R~S,A  (37) 

and  hence 

na  - ~ R3'2.  (38) 

and 

B2  ~ P - Rs/2-  (39) 

In  this  limit  the  beam/plasma  system  decompresses  during  implosion,  with  n.  p and  Bn  decreas- 
ing. 

The  actual  result  obtained  by  numerical  solution  of  the  equations  naturally  lies  between 
these  two  extremes.  The  ring  is  always  observed  to  compress,  but  at  a rate  slower  then  that 
given  by  Eqs.  (34-35),  and  the  scaling  is  not  a power  law  in  R.  If  I,,  = 0 initially,  the  behavior 
tends  to  resemble  the  second  model  increasingly  as  turnaround  is  approached.  The  dependence 
of  the  degree  of  field  reversal  on  the  magnitude  of  the  electron  current  induced  during 
compression151  explains  why  attempts  to  derive  a scaling  law  for  this  parameter13  81  do  not 
appear  to  yield  a simple  result.  There  is,  in  fact  no  clear-cut  way  to  predict  the  scaling  without 
specifying  the  geometry  of  the  compression. 


IV.  COLLISIONS 


The  electron  thermal  spread  is  assumed  to  be  much  larger  than  the  thermal  spread  of 
either  ion  distribution  or  the  relative  drift  between  any  two  species.  The  average  momentum 
transfer  rate  resulting  from  a collision  between  particles  of  species  a and  /3  is  given  by 


where 


m, 


dv„ 


dt 


~ v?/,3ma(v„  - v„) 


(40) 


4nZgZre*  (1  + mB/mT) In  A nT 


mi  vbt 
11 


(41) 
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tyJTn  Zle*(\  + mJm<)mJnnl\n  A 


1 3 ml(kTyn 

a = B,  T,  and,  from  conservation  of  momentum, 


(42) 


r»/0  . 


',awa*'s  (43) 

Here  In  A is  the  form  of  the  usual  Coulomb  logarithm  appropriate  to  the  species  pair  o,/3,  and 


v up 


v„  - vfl 


Correspondingly,  the  average  temperature  rate  of  change  resulting  from  a col- 


lision is. 


dT  r 


dr 

dTT 


dr 


8n  ZgZ^e*  nj\n\ 

t 3 mB  VB7 

_ zBZre 4 ”alnA 

b 3 mT  VBt 


(44) 

(45) 


for  ion-ion  encounters,  and 
dT,. 


dr 


8V2tr  Zle*y[me  n,, InA 


ma(kTe) 


3/2 


k ( T,  - TJ. 


(46) 


a = B,  T,  for  ion-electron  encounters,  with  the  remaining  rates  ( dTJdr )„  defined  so  as  to  satisfy 
conservation  of  energy. 

Consideration  of  the  magnitudes  of  these  rate  formulas  reveals  the  following  general 
features:  (/)  both  electron  and  target  ions  contribute  significantly  to  the  rate  at  which  beam  ions 
slow  down;  (//)  the  relative  velocity  with  which  beam  ions  move  with  respect  to  the  target  ions 
is  chiefly  affected  by  B—T  collisions,  because  electron  collisions  act  in  the  same  sense  (as  a 
drag)  on  both  ion  species;  ( in)  thermalization  of  the  beam  also  results  principally  from  colli- 
sions with  target  ions. 

On  the  basis  of  these  generalizations,  we  can  estimate  the  relative  slowing  down  of  beam 
and  target  ions  through  collisions  as 


— (vB  - vT)C0„  = - (vf/r~  v[/B)  (v8  - vT) 
= — Vs(v  g — v J-)  . 


(47) 
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For  the  usual  case  where  the  target  ion  mass  density  substantially  exceeds  that  of  the  beam. 
nrmr»  nBmB,' Eq.  (47)  implies 

(48) 

at  the  same  time,  the  adiabatic  compression  produced  by  the  imploding  liner  tends  to  cause 
both  ion  species  to  accelerate  in  the  azimuthal  direction  according  to 


dv„ 


dt 


adiab 


V"  R V“  R, 


Taking  the  difference  between  the  beam  and  target  equation  (49)  yields 


d . . R/ 

V B ~ v f)  Adiab  — — ^ (vfl  — V jr)  . 


Eqs.  (47  and  (50)  give  for  the  net  time  rate  of  change  of  the  relative  velocity 


(49) 


(50) 


(jl  ~ vr)  — — ( R//Ri  + r',)  (vfi  - \T) 
The  condition  that  this  relative  velocity  be  a constant  is  thus 


(51) 


*//*/  (52) 

When  Eq.  (57)  is  satisfied,  the  beam  is  said  to  be  clampedno].  With  a tritium  target  there  is  an 

advantage  in  clamping  the  beam  at  a relative  energy  « - ~ mB  v|r  ~ 150  keV  which  maxim- 
izes the  reaction  rate  for  D-T  fusion. 

Clamping  is  of  course  accompanied  by  a monotonic  increase  in  thermal  energy  according 
to  Eqs.  (44-45).  The  ion  thermal  energy  density  w‘,h  = y k(n„TB  + nTTT)  increases  as  a 
result  of  ion-ion  collisions  at  a rate 


dwjh  _ 4 irZlZfe*  nBnT InA 

dt  V gj- 

Using  (48),  we  see  by  comparison  of  (52)  and  (53)  that  the  time  scale  for  implosions  of  the 
liner  is  comparable  to  that  for  heating  up  the  ion  beams.  The  electron  heating  rate  can  be  even 
faster. 


mB  mT 


(53) 
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Note  thai  if  were  approximately  constant,  the  clamping  condition  (52)  would  imply  an 
exponential  decrease  in  with  time.  As  this  is  not  realizable,  clamping  evidently  cannot  be 
maintained  close  to  turnaround. 


In  differencing  the  equations  in  the  code, we  found  it  convenient  to  use  as  dependent  vari- 
ables quantities  that  are  approximately  conserved.  Thus  instead  of  T„  we  used  the  entropy 
functions  lEqs.  (19-21)),  which  now  satisfy  equations  of  the  form 


4 £^(7^-  Ta),  (54) 

v 

where  the  v'f  " are  defined  as  the  rates  in  Eqs.  (44-46).  Similarly,  the  slowing-down  rates  enter 
as 

4 (*vu)  = R Z^Vp-  v„).  (55) 

Ul  p 

V.  OTHER  DISSIPATIVE  PROCESSES 

Collisions,  discussed  in  section  IV,  can  transform  directed  energy  into  thermal  energy. 
Although  essential  for  clamping,  they  may  be  deleterious  if  they  (i)  increase  the  ratio  of  beam 
ion  gyroradius  to  ring  thickness  ettessively;  (ii)  cause  too  much  of  the  liner  energy  to  go  into 
pumping  up  the  target  plasma;  or  (iii)  lead  to  premature  loss  of  confinement  as  a result  of 
decrease  of  beam  current  below  that  needed  for  field  reversal.  In  addition,  the  following  loss 
processes  can  remove  energy  from  the  system  entirely:  radiation,  heat  conduction  along  field 
lines,  particle  diffusion  across  lines,  charge  exchange  with  impurities,  and  ohmic  heating  within 
the  liner.  The  last  of  these  can  have  a second,  more  serious  consequence:  finite  resistivity 
gives  rise  to  diffusion  of  field  lines  through  the  liner,  untrapping  the  magnetic  flux  which  holds 
the  ring  at  a safe  distance  from  the  liner. 

Radiation  processes  are  modeled  by  adding  loss  terms  to  the  expression  (55)  for  the  time 
rate  of  change  of  the  electron  entropy  function.  For  bremsstrahlung  and  synchrotron  (cyclo- 
tron) radiation  we  have  the  terms 
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■Jr  I*"-1  r,)^-  - V*~l  x 5.35  x 10  24  WD  + NtZ})  T} 


■J  Wy-'T')cyi-  - yy-'  X 3.98  X lo16  J~J.  (57) 

where  T is  given  in  eV,  02  - y kTe/m,.c 2 and  B„  - nulr/2irr.  In  the  spirit  of  the  circuit- 

theoretical  approach  (wherein  the  ring  is  a macroscopic  circuit  element  with  certain  lumped 
parameters  derived  from  microscopic  processes),  the  radiation  rates  are  calculated  by  averaging 
the  field  strength  over  the  ring  cross-section. 

In  the  same  fashion,  thermal  conduction  losses  can  be  treated  by  writing 

J,  \yy~'  r“Lw  = " yY~'  4*2Rr  (W-  - 4rr/?  K„Ta.  (58) 

where  kq  is  the  average  cross  field  thermal  conduction  of  species  a.  The  fastest  thermal  loss 

process  is  that  associated  with  the  target  ions,  a = 7.  Furthermore,  thermal  equilibration, 
alpha-particle  heating,  etc.,  can  be  included  in  an  average  sense  in  the  same  form. 

Finally,  particle  losses  can  be  estimated  simply  by  assuming  smeared-out  density  profiles 
according  to  some  law  like  the  Bennett  pinch.  If  a given  profile  extends  past  the  position  of  the 
separatrix,  located  at  average  minor  radius  r = rs,  that  portion  of  the  particles  located  at  r > r, 
is  lost.  A simple  calculation  then  gives  the  loss  rate  as  the  rate  at  which  particles  "fall  over  the 
edge."  Thus  we  find 


dNa  Gr 2 

hf  = - (59> 

l dt  Uff  rs 

where  G is  a geometrical  factor  (equal  to  12  for  a Bennett  profile)  which  decreases  as  the 
assumed  profile  becomes  more  localized,  and  is  the  total  scattering  rate  for  species  a. 


VI.  A NUMERICAL  EXAMPLE 

Using  the  equations  and  techniques  described  in  Section  II- V,  we  consider  the  following 
situation.  A liquid  lithium  liner  (density  p = 0.54  g/cm3)  of  length  L = 13.5  cm  and  inner  and 


BOOK,  TURCH I,  AND  SH  IN 


outer  radii  31.59  cm  and  48.43  cm  implodes  with  velocity  3 x 104  cm/s  on  a fully  ionized  D- 
HeJ  ring  with  major  and  minor  radii  of  30  cm  and  0.758  cm,  respectively.  The  initial  target  ion 
number  densities  are  nHfi  = 2nI}  = 3.59  x 10l6cm  3.  The  temperatures  are  TL j = 23.7  keV, 
T j - 1 keV  and  T(,  = 10  keV.  The  deuterium  current  is  1.52  MA,  twice  the  electron  back 
current.  These  numbers  are  chosen  to  give  a beam  ion  streaming  energy  of  550  keV  and  a 
poloidal  field  of  200  kG,  with  beam  clamping.  Since  the  emphasis  was  on  determining  £),  only 
the  part  of  the  evolution  in  the  vicinity  of  liner  turnaround  was  considered,  and  the  early-time 
conditions  giving  rise  to  these  parameters  were  not  investigated. 

Figure  3 shows  how  the  beam  and  liner  radii  change  in  time.  Note  that  the  separation 
increases,  a reflection  of  the  increase  in  beam  minor  radius  (Fig.  4).  Correspondingly  the 
number  densities  (Fig.  5)  drop,  level  off  as  collisional  heating  and  compression  come  into  bal- 
ance, then  drop  again  in  the  expansion  (decompression)  stage,  and  the  poloidal  field  (Fig.  6) 
decreases,  increases,  then  decreases  monotonically  after  turnaround.  The  various  forms  of 

energy  (magnetic,  liner  kinetic,  ion  directed,  and  thermal)  are  plotted  in  Fig.  7,  along  with  the 
fusion  yield.  Figure  8 shows  how  the  component  temperatures  increase  near  turnaround,  the 

evident  irreversibility  being  a consequence  of  collisions. 

Running  time  on  the  calculation  using  an  IBM  360/168  was  91  seconds,  of  which  about  a 
quarter  was  required  for  diagnostics.  Using  ten  current  loops  to  represent  the  liner  current 
profile  approximately  doubles  the  running  time,  since  roughly  twice  as  many  differential  equa- 
tions have  to  be  solved  [matrix  inversion  of  Eq.  (2)  does  not  add  any  substantial  amount  to  the 
total].  It  turns  out  to  be  convenient  in  writing  the  code  to  make  extensive  use  of  nested 
sequences  of  statement  functions  in  redetermining  force  balance  on  each  time  step,  and  most  of 
the  running  time  is  expended  in  this  task. 

A variety  of  prescriptions  are  possible  for  defining  the  initial  conditions.  The  main  thing 
is  to  insure  that  they  be  neither  over-  nor  underdetermined.  When  working  with  multiple  liner 
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current  loops,  we  arbitrarily  imposed  the  condition  that  the  flux  threading  all  the  loops  be  the 
same.  Though  straightforward,  this  is  unlikely  to  be  a good  approximation  in  the  late  stages  of 
the  implosion  if  finite  liner  resistivity  is  modeled. 

VII.  CONCLUSIONS 

We  have  presented  a new  numerical  technique  for  solving  problems  involving  the  dynam- 
ics of  charged  particle  rings.  Its  principle  advantage  is  that  it  is  couched  in  circuit-theoretical 
terms,  obviating  the  need  for  solution  of  partial  differential  equations.  Because  of  its  adaptation 
to  the  physics  and  geometry  of  such  problems,  the  method  can  be  implemented  with  only  a 
small  number  10)  of  current  carrying  elements.  In  effect,  it  replaces  the  uniform  or  quasi- 
uniform mesh  of  the  standard  2D  finite-difference  technique  with  a highly  nonuniform  "mesh" 
of  circuit  elements,  located  optimally  to  reflect  the  relevant  physics. 

The  code  has  been  applied  to  calculations  of  the  thermonuclear  yield  and  other  charac- 
teristics of  a beam-target  fusion  device.  The  particular  concept  for  which  the  code  was  originally 
developed  turns  out  to  be  disappointing  in  terms  of  its  efficiency  as  a reactor  (the  examples  of 
Section  VI  yielded  Q ~ 3.2%),  and  also  appears  to  be  unstable  to  kink  modes";  however  it 
may  have  non-fusion  applications.  It  is  clear  that  the  code  can  be  applied  to  a variety  of 
axisymmetric  situations  involving  field  reversal  and  changes  of  system  geometry,  and  therefore 
is  potentially  of  wider  utility. 

This  work  was  supported  by  the  Office  of  Naval  Research. 

References 

1.  P.  Dreike,  R.  L.  Ferch,  A.  Friedman,  S.  Humphries,  R.  V.  Lovelace,  G.  Ludwig,  R N. 
Sudan,  D.  L.  Book,  G.  Cooperstein,  A.  T.  Drobot,  J.  Golden,  S.  A.  Goldstein,  C.  A. 
Kapeianakos,  R.  Lee,  W.  Manheimer,  S.  J.  Marsh,  D.  Mosher,  E.  Ott,  S.  S.  Stephanakis, 


kL 


17 


BOOK,  TURCHl,  AND  STEIN 


and  P.  J.  Turchi,  Plasma  Physics  and  Controlled  Nuclear  Fusion  Research  1976  (IAEA, 
Vienna,  1977),  v.  3,  p.  415. 

2.  S.  J.  Marsh,  A.  T.  Drobot,  J.  Golden  and  C.  A.  Kapetanakos,  Phys.  Res.  Lett.  39,  705 
(1977). 

3.  R.  N.  Sudan  and  E.  Ott,  Phys.  Rev.  Lett.  33,  355  (1974). 

4.  D.  E.  Baldwin  and  T.  K.  Fowler,  Lawrence  Livermore  Laboratory  Rept.  UCID-17691 
(1977). 

5.  D.  E.  Baldwin  and  M.  E.  Rensink,  Lawrence  Livermore  Laboratory  Rept.  No.  (submitted 
to  Phys.  Rev.  Lett.). 

6.  T.  Ohkawa,  Nuclear  Fusion  10,  185  (1970). 

7.  D.  L.  Book  and  P.  J.  Turchi,  Bull.  APS  20,  1278  (1975).  Abstract 

8.  E.  S.  Weibel,  Phys.  Fluids  20,  1195  (1977). 

9.  S.  I.  Braginskii,  Reviews  of  Plasma  Physics , vol.  1 (Consultants  Bureau,  New  York,  1965), 
p.  205. 

10.  10.  J.  M.  Dawson,  H.  P.  Furth  and  F.  H.  Tenney,  Phys.  Rev.  Lett.  26,  1156  (1971). 

11.  R.  N.  Sudan  and  M.  N.  Rosenbluth,  Phys.  Rev.  Lett.  36,  972  (1976);  R.  N.  Sudan  and  M. 
N.  Rosenbluth  (to  be  published). 


NKL  MEMORANDUM  REPORT  3827 


e 

C ’ION  8EAM-PLA$Ma  INTERACTION  CLAMPED  THROUGH  AXIAL  C OPPRESS  I ON" 

c 

C this  program  calculates  compression,  heating  and  burn  in  an  axially 

c compressed  beam. target  plasma,  the  beam  is  a ring  or  o ions,  the 

C TARGET  a MIXTURE  OF  COLO  TRITIUM  AND  HOT  ELECTRONS.  the  CONFINING 

C MAGNETIC  FIELD  BUILDS  up  OPING  to  COMPRESSION  produced  by  an 

C IMPLODING  LITHIUM  LINER  (LINUS),  DRIVEN  ELCTROMAGNETicALLY  BY 

C external  theta-pinch  findings,  parameters  are  chosen  so  that 

c COMPRESSION,  C8LLISI0NAL  SLOPING  OF  THE  D BEAM,  ANO  LOSSES  BALANCE, 

c producing  a situation  in  phich  the  beam  is  clamped  for  the  entire 

C OURATION  OF  THE  pOMPRESSION, 

C 0.  L.  BOOK,  P.  J.  TURCMI  A D.  L.  STEIN 

c 

IMPLICIT  REAL»8  (A-I,  K-Z),  INTEGER**  (J),  REAL**  (S) 

LOGICAL**  COLL.  LOSS,  BURN,  FORPRD 
COMMON/  NEPTON/  DTO,JITER 
common  /G8L0CK/  2(25) 

COMMON/  array/  Y (2S ) , DY (25) , Y0(25) 

EXTERNAL  DERI V 

NAMELIST  /liner/  Rio,  R20,  RP,  B20,  *Hfl,  OMEGA,  LNGTHO,  ETA, 

1 SPHT,  TL,  PCA 

DATA  RIO,  R20,  RW,  B20,  RHO,  OMEGA,  LNGTHO,  ETA,  SPHT,  TL,  PCA 

1 / 2.0602,  ?. 12602,  J.D2,  O.DO,  .5«Do,  0.00,  Ij'.500, 

2 *5 . 2D-7,  6.07,  5.D2,  .IDO  / 

NAMELIST  /PLASMA/  RMAJOR,  Pd,  TO,  TE,  TT,  BP,  VRATIO,  NRATI8, 

1 1RATI0 

DATA  RMAJOR,  WD,  TO,  TE,  TT,  BP,  VRATIO,  NRaTIO,  IRaTIO  / 

I 2. 02,  7.D2,  i .02,  I.D2,  1.02,  S.D5,  .500,  I. Do,  .500  / 

NA“ELIST  /LOGICL/  COLL,  LOSS,  BURN,  FORPRD 
OATa  COLL,  LOSS,  BURN,  FORPRD  / .TRUE.,  '.TRUE.,  .TRUE,, 
l .FALSE.  ( 

NAMELIST  /CNTRPL/  DT,  DTPLST,  DTFILM,  DTDUMP,  TLAST,  T,  TPLOT, 

1 . TFILM,  TOUfP 

DATA  OT,  DTPLOT,  OTFILM,  OTDUMP,  TLAST,  T,  TPL8T,  TFILM,  TDUMP 
) /l.D-5,  l.O-a,  1,0-),  1.0-2,  .100,  0.00,  -1.0-12,  -1.0-12, 

1 -1.0-12  / 

CALL  INDUMP 
CALL  CRDUmP (40,10) 

e 

C INITIALIZE. 

c 

to  READ  (5,LINER,END»t00) 

READ  (5, PLASMA, ENDpIOO) 

READ  (5, LOGJCL, END»100) 

RE*0  (5,CNTR0L.EN0«100> 

PRINT  *B 
DO  11  J * 1,  5 

tl  PRINT  RO 

PRITE  (6, LINER) 

WRITE  (6, PLASMA) 

PRITE  (6, LOGICL) 

PRITE  (6.CNTR0L) 

00  IS  J « 1,  25 
15  0(J)  ■ O.DO 

JITER»0 
DTMAX  ■ OT 
DTBpDT 

CALL  INPUT (T, 

1 RIO,  R20.  »P,  B20,  RHO,  OMEGA,  LNGTHO,  ETA,  SPHT , TL,  PCA, 

2 RMAJOR,  PD,  TD,  TE,  TT,  BP,  VRATIO,  NRATIO,  IRATIO, 

) COLL,  LOSS,  BURN,  FORPRD, 

T P n % 1 T r OTPLOT,  OTFILM.  OTDUMP,  TL»»T,  TPlBT,  TFILM,  TDUMP) 
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C 

e PERFORM  0UMf*8  * DIAGNOSTICS  AT  APPROPRIATE  INTERVALS  * 

e 

20  IE  (T  .lT.  TDUHP)  GO  TO  JO 

TDUHP  * TDUMP  ♦ DTDUMP 
JO  IE  (T  .IT.  TEILM)  GO  TO  40 

Tp Jl*  ■ TFtlM  t DTEILH 

40  IE  (T  .IT.  TPlOT)  GO  TO  50 
TPL8T  • TPLOT  ♦ OTPLOT 
CALL  DERIVfT, JTOTAL) 

IE  (T  ,E0.  0,00)  PRINT  40 
IE  CT  ,NE . 0.00)  PRINT  41,  T 

OATA  ST,  JYO.  SOY,  SLI,  SRINO,  |PLAS.  INU,  ICNU,  SOATA  / 

1 'V  'TO  'DY  ',  'Ll  ',  'RING ' , <PLA8»,  *NU 

•I  'CU  'DATA'  / 

CALL  OSPLAY(SY) 

CALL  DSPLAYCSTO) 

CALL  OSPLAY(SDY) 

CALL  DSPLAY(SLI) 

CALL  DSPLAY ( SRING) 

CALL  DSPLAY(SPLAS) 

CALL  DSPLAY ( INU ) 

IE  ( T .GT.  6.00)  GO  TO  50 
CALL  DSPLAY ( ICNU) 

CALL  OSPLAY(SDATA) 

SO  CONTINUE 

c 

C ADVANCE  VARIABLES  ONE  TIMESTEP. 

c 

JITERro 

IE  (Y(J1  ,N£.  6.00)  DT  ■ DMINKOTNAX,  . 6500*DABS(Y(2)/Y(3)}) 
DTOrDT 

CALL  INT(JT0TALiT,OERIV,DT) 

IE  (T  .LT.  TLAST)  GO  TO  20 
GO  TO  10 

oo  pormatco',  'initial  display,  including  tabulated  constants!') 

41  FORMATCM',  ' 0 1 SPLAY  AT  T JMg  T R 1PDI0.4,  ' SECl') 

46  ESRMATl't',  to*,  M M t A C //  ' 3SX , *IBN  BEAN-', 

l 'PLASMA  INTERACTION  CLAMPED  THROUGH  AXIAL  COMPRESSION') 

44  FORMAT!'*',  J5X,  'I  ', 

1 'P  I C AC  ') 

100  RETURN 

END 

C 

SUBROUTINE  DERIVUJME, JTOTAL) 

IMPLICIT  REAL«8  (A-I,  K-Z),  INTEGERS  (J),  REALaR  (S) 

INTEGER*®  MOO 

DIMENSION  Rl INER (5) » EMAGK5) 

OIMENSION  ZA(J,J),  ZR(J) 

DIMENSION  Rmaji (#),RMAJ2(4),RMIN1 (4),RMlN2(fl) 

COMMON  / NEWTON  /DT,  JITER 
COMMON  / JNCO  / MU,  KMU 
COMMON/  RAOII/  RMInOR,RMAJOR,RJ 
COMMON  / PBLOCK  / P 
COMMON/  COEFE/  KCE.NDTOT,  RVD 
COMMON  / CURREN  / ILINERCJ) 

COMMON/  LENGTH/  LnGTh  C5 ) , LHLSQC5) 

COMMON  / INDEX  / JMAX,  JMAX1 

common/  flux/  phi,psii(5) 

COMMON/  ARRAY/  YC25),  0Y(25),  Y0(25) 

LOGICAL*®  NOCOLLi  NOLOSS,  noburn, FORMRD 
OIMENSION  SLINE(JJ) 

c 

c the  following  statement  punctions  are  used  eelow  in  solving  for 

C RMAJOR  and  RMINoR  by  RfGUIRING  THAT  THE  FORCES  IN  THE  CORRESPONDING 
C OI«ECTIONs  VANISH.  MERE  A,  B,  C A®E  RESPLCTIVELY  THE  MJNOR  *ND 

C MAJOR  RADII  OF  THE  RING  AND  THE  LINER  INSIDE  RADIUS,  ALL  IN  CM. 
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c 

LR(A.B)  « KHU«B»(OL8G(6.00»B/A)  - 2. DO  ♦ KD) 

D *LR  t * « B ) ■ -KMU*8/A 

DBLR(A»B)  ■ KMU*(0L0G(fl,D0»8/A)  - 1.C0  ♦ KD) 
LL(C,D)*HU*C*(DL0G(16.DO*C/D)  - 0.5D0) 
DCLL(C,D)«HU*(0L8G(lt.D0»C/0)  ♦ o.SDO) 

LN(D,F)»0.5d6»  (D/F)»DLOG(D) 

P AC (0,F)«(((F-0)*»2)/C 2.00*0 «F))*DL8G(DABS(F-D)) 

Hll (e,0,F)*HU.C*(DLBG(16.D0*C)-LN(O,F).LN(F,D)*FAC(D.F)-0.SD0) 
DCHLL(C,D,F)«HU»(DL8G(16.D0»C)-LN(0,F).LN(F,0)*F*C<D.n*0.500) 
HLRC8fC,F)»MU*DSQRT(B»C)*(.5D0*DL8G(6.4CI*B*C/(CB-C)«»2  ♦ 

1 F))  . I. 00  - ((B-C)/DS0RT(F)),DATAN(DS0RT(F)/(8-C))) 

DB'n.R(B,C,F)*.  JD0»MU*DS3RT  (C/B ) « ( .500 *DL8G (6 . 40 1 *B *C/ ( CB-C ) *«2 
1 ♦F))*((C-J.D0*B)/DS0RT(F))*DATAN(D8GRT(F)/(B-C))) 

OCHLR(B,C,F)*.5DO»MU«OSQRTCB/C)*{.5DO»OL8G(6.401«B.C/( (B-C)**2 
1 ♦F))*(Cj'.DO«C-P)/DSQRTCF))*0*TAN(DSCRT(F)/CB-C))) 

HtR*f*,0.C,O,F)»MLR(B,C»F)*(OSORT(LR(A.B)*LL(C.D))-HLR(B,C,F)) 
1 /(I.D0*((C-B)/A)**P) 

D*HtRj(A,B,C,D.F)«(.5D0»DALR(A,B)»D8QRT(LttC.0)/LR(*.B))*(R/*) 

1 «((C«B)/A)**P*(DSQRT (LR(A,8)*LL(C#D))-MLR(8,C,F))/(1 .00* 

2 ((C*8)/A)**P))/(1.00  ♦ ((C-B)/A)»»P) 

DBHLR$(A,B.C.D.n»nBHUR<8,C,F)A{.5D0«DBUR(**B)*D8ORT(l.L(e»D)/ 

1 LR(A.B))  ♦ (P/(C-e))*((C-0)/A).*P*(O3ORTCLRCA,B)*LL(C,O))- 

2 HLR(B,C,P))/'(l.DO  ♦ ((C*B)/A)**P))/(1 .DO  ♦ ( (C-B ) /A ) **P ) 
OCHI.R$CA,B.C.D.F)»DCHLR(B,C,F)*(.500*DCLL  (C , 0 ) »0SQRT  (LR  (A,  6)  / 

1 LL(C.D))-(P/jtC-B))*CCC-B)/A)*.P*(DSQRT(l.R(A,B  )»LL(C,D))- 

2 hi.R(B,C,F))/(1.D0  ♦ (CC-B)/A)**P))/(1,D0  ♦ ( (C«B)/A)*«P) 
NKT(A.B)  ■ PVG/(KV8L*A*A*B)««GHl 

8 IGMA (E ) ■ l.D-24*(A5  ♦ *2/(1. Do  * (AJ*E  - A«)**2))/(E* 
l (0EXP(*l/080RT(E))  - 1 .00) > 

DATA  ET8T0  / 0.00  / 

c 

C TIME  (INDEPENDENT  VARIABLE)  HA8  UNIT  DERIVATIVE. 

C 

D V f 1 ) ■ 1.00 

c 

C REPLACE  SUBSCRIPTED  QUANTITIES  WITH  MORE  FAPILIAR  NOTATION, 

C 

RtSQ  ■ Y (2) 

RU  * Y C 3 > 

PS  12  « Y (5) 

RVD  ■ Y (7) 

RVT  » Y (8 ) 

TDVGH 1 ■ Y (9) 

TEVCM 1 ■ Y ( 1 0 ) 

TTVGHJ  • Y ( 1 1 ) 

NDTOT  ■ Y ( 1 2 ) 

NTTOT  ■ Y(]J) 

E8HMIC  « Y(lfl) 

CRAO  • Y ( 15) 

NC8UNT  ■ Y ( It ) 

PVOPRD  ■ Y ( 1 7) 

BURNUP  ■ Y ( 1 6 ) 

D8  5 Jil.jHAxi 
5 PSU  (J)«Y(19tJ) 

PHl«Y(19*JHAX) 

c 

C PIND  88HE  OF  THE  QUANTITIES  NEEDED  TO  DESCRIBE  THE  LINER  DYNAMICS, 
C 

psouso  ■ ru*ru 

R28Q  > R 1 SO  * P 080 
PI  ■ DSOPT(RISO) 

R2  ■ DSQFT (R2SQ) 

U1  » RU/R 1 

n I nfrPT  U2  « RU/P2 
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C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 
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~ ' 8l  a KFLUX*PSli  (D/R18Q 

82  a KFLUX*PSI2/(RW30  - 8280) 

RISC  a DLOG (R2SQ/R1 3Q) 

DRSO  a R 1 OSS  • RISC 

P12SQ  ■ I.Do/RtSO  - I.D0/R2SQ 

PMAG1  ■ XP*61*B1 

PMAG2  • KP«B2*B2 

U1BYR1  a Ut/Rl 

NO*  DO  ION  RING  DYNAMICS  AND  PLASMA  PROCESSES.  USE  ALL  THE 

available  algebraic  relations  before  computing  any  derivatives, 
the  ring  major  and  minor  radii  are  founo  using  newton's  method  to 

SOLVE  the  EQUATIONS  FOR  FORCE  BALANCE  IN  THE  RING*  GIVEN  THE  PLUXe» 
WHICM  AHE  ENCLOSED  BY  THE  RING  and  liner  (Phi  ano  psii, 
RESPECTIVELY),  AND  USING  HANDBOOK  FORMULAS  FffR  SELF-  AND  MUTUAL 
INDUCTANCE*.  MK8  UNITS  ARE  USED  IN  THIS  PORTION  OF  THE  COOE. 

EXTRAPOLATE  FROM  LAST  TWO  TIME  STEPS  TO  GET  GOOD  INTIIAL  GUESSES 
FOR  RMAJOR,  RM I NOR  (USEO  ONLY  ON  EVEN  STEPS  OF  R-K-G). 


IF  (JITER  .NE.  0)  GO  TO  1 

JSTEP  • 0 

0T3L0  « DTN£W 

OTNEH  « DT 

NT!  » -DTNEX/DTOLD 

KT2  ■ 1.00  - KTt 

CONTINUE 


JSTEP  » JSTEP  ♦ 1 

IF  (MOD {JSTEP, 2)  .NE.  0)  GO  TO  2 

RMAJOR  ■ WT1*RNAJI (JSTEP)  ♦ NT2*RMAJ2CJ3TEPJ 

RMINOR  ■ WT  J 0RMIN1 (JSTEP ) ♦ WT2»RMIN2 (JSTEP) 

CONTINUE 


JIT  « 0 

A ( RMINOR 

B ■ RMAJOR 

RJ  m R 1 

PvG  * KNKT*(NDTBT»(TDVGM1  * ZD*TEVGMI)  ♦ NTTOT* (TTVGM l a 
I ZT*TEVGM1))  t 1 ,D-7 »PvGPRD 

CONTINUE 

DAFA«(PA(A*0/2'.D0,B,RJ,NKT(Aa0/2.0O,8));fA  ( A-0/2 .00  , b ,R  J , 

1 NKT(A-0/2.p0,B)))/0  , . 

DAPB«(FB(A*0/2.DO,B,RJ,HKT (AAO/2.DO,B))iPB  < A-0/2 , DO , B , R J , 

1 NKT(A.0/2.D0,8J))/0 

OBPA»(FA(A,b*0/2.Do,RJ,NKT(A,Ba0/2.D0)).PA(A,B-0/2.OO,RJ, 

l N*T(A,B-0/2.D0)))/0  . 

DBP8»(FB(A,8a8/2.DO,RJ,NKT(A,BaO/2.00))-P8(A,B-0/2.00,RJ, 

I NKT(a, 8-0/2. 00))  )/8 

DET  * 0APA»DBF8  - DAFO-DBFA 
FMINOR  a FA(A,8»RJ.NKT(A,8)) 

FMAJOR  ■ F8(A,8,RJ,NKT(A,B)) 

DA  » (DBP8»FMJN8R  - DBF  A *FMa JBR ) /OET 
OB  a (OAF a *PN^ JOR  - DAF8*FMIN8R)/qeT 
A ■ A • DA 

8 ■ B - OB 

JIT*JIT»l 

IF  (0AB3CDB)  , 5T.  RTESTaB  ,8R.  DABS(DA)  .GT.  RTEST*A)  00  TO  10 

RMIN8R  * A 

RMAJ8R  a 8 

JITER»JIT£RtJIT 

ITER  » DFLOATCJITER) 

IF  (MOO (JSTeP, I)  .NE.  0)  GO  TO  12 
RMAJt (JST£P)  a RMAJ2(JSTEP) 

RMAJ2 (JSTEP)  a RMAJOR 
RMIM(JSTEP)  a RMIN2(JSTEP) 

RmIN2(JSTEP)  a RMINSR 
CONTINUE 
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DEFINE  THE  F9U8WINQ  FAR  DIAGNOSTIC  PURPOSES’. 


umiCRJ.LNCTHCtn 

LL2»LL<RJ*LNGTH(2)) 

MLJ2*HU(RJ,LNGTH(2),LNGTH(n) 

MLRSI»HLRSCRNInOB,RHAJOR,RJ#LNGTH(U,LHL8QCI>) 

EVERYTHING  ELSE  CAN  NOW  BE  CALCULATED. 

LRING  » LR (RH IN8F , RHAJBR) 

LLlNER«ll(Pj,LNGTH(JNAXl)) 

IRINGiILINER(JHAX) 

VO  • RVD/RHA J8R 
VT  » RVT/RHAJ8B 

V8LUH£  ■ KV8L*RHAJ8R*RHIN8R»«2 

NO  « NDT8T/V8LUME 

NT  ■ NT T8T / VPLUME 

NETST  a NDT8TA70  t NTT8T*2T 

NE  ■ NET8T/V8LUME 

10  « KI0«NDT8T*V0/RHAJ8R 

JT  ■ KIT*NTT8T»VT/RHAJ8R 

IE  • IRING  • 10  - IT 

VE  « RMAj8R«It/ £KIE«NET8T) 

HAVING  C8HRUTE0  TMr  magnETBOYNAHU  PARAMETERS#  HE  CAN  WRITE  DOWN 
The  eguatisns  bf  mbtion  8F  The  liner. 

CYC?)  ■ 2 , 00  »RU 
JH*x2ijy«xt.i 

F8RCE»0,00 
08  14  J=t,JHAXl 

F8RCE»F8RCE»ILINER(J)*(0.500*ILINER(J)«CCLL(RJ,LNGTh(J)) 
l ♦IRING»OCHlRS(RMIN8R,RMAJ8R,RJ,LNGtw(J),LHLaO(Jin 

IF(J.EQ.JMAXl)  G8  T8  14 
JP1»J»1 

08  IS  JJ«JPl,JMAXt 

F8RCE«F8«CE+ILINER(J)*ILINERCJJ)*0CHLL{RJ.LHL8O(JJ)fLHLSO(J)) 

C8NTINUE 

Pi  » KFUFBRCE/RJ 
»2  ■ PMAG2 

DY(S)  ■ CRSGUSObR 1 2S0  ♦ SHEGA**2* (R03G  < DRSO« C2 .00*RL8G  * 

1 0RS0«R123G)J  ♦ 2.00*(Pl  • P2)/RH8}/RL8G 

PINO  TEHPERATURES  FR8H  PRODUCT  BP  V T8  POWER  OAMHA  • J AND  T. 

VGh 1 s V8LUHE  **GH 1 
TD  « TDVGMl/VGMl 
TE  » TEVGH 1 / V C M I 
TT  * TTVGH1/VGH1 

KEEP  TRACK  8T  RING,  PLASHA  ENERGETICS. 

EKIN  a Krf  IN.R35US3.RLBG 

erst  ■ kpbt*(orso**2*rlbg/2,do  ♦ roso»cdrsg  a criso  ♦ r2sgj/ 

l 4.00)) 

EHAG2»LNGTH(jMAXn •B2»62*(RW80-R250)/8.D0 

EHAG«0.507*IR!'.G«PHl 

08  lfc  J«1,JWAX1 

EMAGI IJ)*0.5O7»ILINER(J)»P8ll (J) 

EMAG*EHAGaEHAG1 CJ) 

CONTINUE 
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EMAG»EHAG*EMAG2 

I ynt/81  "0  " HP*V0*»?/(2.00*B0LTZ) 

»T  ■ MT»VT»«2/C2.D0»B0LTZ) 

EOOIR  ■ NDT0T*B0LTZ«HD 

ETDIR  ■ NTT0T*B8LTZ*NT 

EDTH  ■ KE*NCT6T*TD 

ETTH  ■ KE*NTTOT*TT 

EETH  » KE»NET8T*TE 

EBEAM  c EDOIR  ♦ EOTH 

EPLAS  « E TP  IP  ♦ ETTH  ♦ EETH 

ERING  « EBEAH  ♦ EPLAS 

EP8T»EBJNC*£HAG 

ET8T  « EKIN  ♦ ER8T  ♦ EP8T 

IE  (ET8T0  .EC.  O.OO)  ET8T0  « ET8T 

YIELD  ■ KYIELD«BURNUP 

0 ■ VIELD/ET8T0 

ENET  ■ ET8T  ♦ ERRD  ♦ ESHHIC  • YIELD 

PD  ■ B8LTZ*N0»T0 

PE  ■ BOLTZ«NE«TE 

PT  ■ B8LTZ»NT*TT 

PT9T  ■ PO  ♦ PE  ♦ PT 

BETA  « KBETA*DSQRTfTE) 

BP  « ,2D0*IRING/'PHIN8R 
KAPPA  ■ KPR8P«VD/(BHIN8R«BP) 

PP8L  » KP«BP*«2 

BET APL  • NPT(R«INOR»  HMAJORJ/PPOL 
TL  ■ E8BHIC/MTCAP  ♦ TlO 
C 

C PINO  ACCELERATION  AT  INNER,  OUTER  PACE  OP  LIN£R, 

C 

81  ■ 8MEGA«Rl6SQ/RlSO 
02  ■ 9MEGA*R20SQ/R23Q 
Cl  ■ COY (3 ) • Ul»Ul)/Rt  - 0J»01«R1 
62  * C D Y (J)  • U2*U2)/R2  - 02«02*R2 
C 

c zero  derivatives  se  isn  ring  and  plasma  quantities  conserved  in 
c the  ABSENCE  be  dissipation, 

c 

JT8TAL*25 
08  15  J»«,  JT8TAL 
19  OY(J)  ■ 0 , Do 

JT8TAL«22 
VDT  ■ VD  - VT 
IE  (NOCOLL)  G8  TB  20 
C 

c PUT  IN  C8LLISIBNAL  CPPECTS,  IE  ANY,  start  by  STARING  THE  RELATIVE 
C VELOCITIES  AND  THEIR  SQUARES, 

c 

VOE  • VO  - VE 
VET  ■ VE  - VT 
VTO  ■ -VDT 
VEO  ■ -VDE 
VTE  « -VET 
V0T3Q  ■ VDTuVOT 
V0E3Q  ■ VDE»VDC 
VETSQ  ■ VET«VET 

c 

NOMD  ■ ND *Hq 
NEHE  ■ NE*Hf 
NTHT  ■ NT«NT 
C 

C PINO  COULOMB  LOGARITHMS. 

C 

L8GLDT  ■ LOGOTO  - . 5D0»DLBG (N£/ (TE«VDT3Ca«2) ) 

EL<*G  ■ ,5D0«DLBC(NE/TE»*3) 

L06LOE  ■ LOGOEO  - ELOG 
L8GLTE  ■ LOGTEO  - ELOG 
C 
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C CALCULATE  THE  EfFECTS  9F  COLLISIONS  BETWEEN  BEAR  AND  TARGET  IONS. 
C THe  L8W-TEMPEHATURE  LIMIT  OF  TRUBNIKOV'S  FORMULAS  l<  UBeO, 

c 

NUSOT  a CNUSDT»L8GLDT»NT/DAB8(VDT#a3) 

NUSTO  a NU8pT*NDMD/NTMT 

OTDST  a CNUTDT*L6GLDT»NT/DAB8(VDT) 

NUTDT  a OTDST/TD 
OTTSD  a OTO$T*NOMD/NTMT 
NUTTO  a DTTSD/TT 

OW08T  a OTDST  - KTD*NU30T*VD*VDT 
OWTSO  a OTTSD  - KTT*NUSTD*VT*VTD 
IF  CU1  .NE.  O.OO)  RNUBYU  a RMAJ6R»NU8DT/U1 
C 

C USE  8LWW-BEAW  LIMIT  TO  COMPUTE  INTERACTION  BETWEEN  I0N8  AND  NOT 
C ELECTRONS, 

C 

EFACTR  a NE/TE««1.500 
NUSOE  a CNUSOE«LOGLDE»EFACTR 
NUTOE  a CNUTDE*LOGLDE*EFACTfl 
OTOSE  a NUTOE«(TE  - TO) 

NU8TE  a CNUSTE*LOGLTE*EFACTR 
NUTTE  a CNUTTE*LOGLTE*EFACTR 
OTTSE  a NUTTEatTE  - TT) 

C 

C CONSERVATION  OF  MOMENTUM  GIVES  THE  REMAINING  NUS'S. 

C 

NUSEO  a NUSOEaNOHQ/NEME 
NUSET  a NUOTEaNTMT/NEHE 
C 

C USING  KNOWN  NUS'S  AND  OT'S,  THE  CORRESPONDING  OH'S  ARE  FOUND. 

C 

DWi;iE  a OTOSE  - KTCNUSDEaVORVOE 
DWTSE  a OTTSE  - KTT*NUSTEaVT*VTE 
C 

C CONSERVATION  OF  ENERGY  GIVES  THE  REMAINING  THREE  OH'S. 

C 

OWESO  a •DWCSEfNO/NE 
OWEST  a -OWTSEFNT/NE 

e 

C FROM  THE  DERIVED  NUS'S  AND  DW'S,  FIND  THE  REMAINING  OT'S. 

C 

OTESD  a OWESO  * KTE»NU8ED*VE«VED 
NUTED  a OTESD/TE 

DTEST  a OWEST  ♦ KTE«NUSET*VE*VET 
NUTET  a OTEST/TE 

C ADO  COLLISIONS  CORRECTIONS  TO  THE  DERIVATIVES  CALCULATED  ABOVE. 

e 

OYtJTOTAL)  bKPHI«RMAJ0R»<NUSED*VEDtNU8ET»VET) 

OY (7)  a -RMajbr»(NU5DE*V0E  ♦ NUSPTaVDT)  - KVOaDY (4 f OTAL ) 

DY(6)  a •RMAJBRatNUSTDaVTO  ♦ NUSTEaVTE)  - KVTaOY  (JTOTAL ) 

DY(R>  a VGM1 • (OTOSE  ♦ OTDST) 

0Y(10)  a VGH 1 • (DTESO  ♦ DTEST) 

DY ( 1 1 ) a VGM1 • {OTTSD  ♦ OTTSE) 

20  IF  (NfiLOSS)  CO  TO  30 

C 

C PUT  IN  OTHER  SOURCES  OR  SINKS.  IF  ANY.  START  WITH  OHMIC  LOSSES. 

C 

OHMICaO.OO 

PCONSTbKR£S»ETA/RLOG 
00  25  Jai.JwMJ 
RLINER(J)«RC0NST/LNOTH(J) 

OHMHTaPLlNER(J)»ILlNER(J) 

OY(lR*J)a*OHMHT 
25  OHMICBOHMICaOHMHT 

C / nJ-PT  DV  (5)  a OPS  12 
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“ 1 ornu)  « o.co 

09  26  J«t,JMAxi 

26  0Y(1«)  » DY£1<()  ♦ t.07»OHMMT«|lIN£6{J) 

c 

c then  treat  lasses  by  diffusion,  charge  exchange,  etc. 
c 

0Y(10)  ■ DY(IO)  ♦ TEL0S3 
DYC12)  « NOLOSS 
0Y<1S)  « NTL0S8 

c 

C FINALLY,  compute  BHEMSSTRAHLUNG  AND  cyclotron  radiation  losses. 

c 

PCYCL  • KCYCL*(BET*»BP)»«2/U.D0  - BETA. *2) 

PB«EM  ■ K0BEm»OSORT(TE)»(ND*ZO»ZO  ♦ NT«ST«ST) 

PRAD  * P9REM  ♦ PCYCL 
DTERAD  ■ KT.PRAD 
NURAO  » OTERAO/TE 
OY ( 1 5 ) » NETOT*PRAD 
DYflO)  s DY ( 1 o ) • DTERAO.VOMI 
30  IF  (NBPURN)  GO  TO  40 

c 

C COMPUTE  THERMONUCLEAR  BURN. 

c 

RATE  * SIGM*(1.D3«WD)*DAB8(VDT)»ND6NTT0T 

0Y(  12)  « OY ( 1 2 ) - RATE 

DY (1 3 ) * OVflJ)  - RATE 

0 Y ( 1 6 ) « KNEUTR»RATE 

DY  1 1 T ) ■ XBURN.VGHI *RATE 

OY ( 1 8 ) » RATE 

60  IF  (1  .GT.  0)  RETURN 

C 

00  35  J»tl,JTOTAL 
35  OY(J)»-OY(J) 

RETURN 

c 

C PASS  in  Data  needed  to  initialise  Run,  saving  initial  values  as  To. 
C PRECOMPUTE  CONSTANTS  FOR  USE  IN  SUBSEQUENT  CALLS  TO  OERIV.  CCS 
C UNITS  ARE  USED  THROUGHOUT,  EXCEPT  THAT  TEMPERATURE?  ANO  PARTICLE 
C ENERGIES  ARE  IN  K£V  AND  CIRCUIT  QUANTITIES  (CURRENTS,  FLUXES  AND 

C INOUCTANCES)  ARE  IN  MKS. 

c 

ENTRY  INPUT(TI“f, 

1 R 1 0 , R20,  RH,  B20,  RMOO,  OMEGAo,  LNGTHO,  ETAO,  SPHT, 

2 TLOO,  pca, 

3 RMAJO,  oDq  , TOO,  TEO,  TTO,  BPO,  VRATIO,  NRATIO , IRATIO, 

6 COIL*  LOS?,  BURN,  FWD, 

3 DTPLOT,  DTFILM,  DTDUMp,  TLA3T,  TPLOT,  TFILH,  TDUMP) 

LOGICAL *4  COLL.  LOSS,  BURNED 

forwrd»fwo 

PHO  » RH90 

OMEGA  * On£ 5 A o 
ETA  ■ ETAO 

TlO  ■ TLOO 

PMiJOR  | RHAJO 
P»l ,Dl 
JH  A X » 3 
JT5T  *L»22 
JMAX1«JMAX«1 

DLNGTH»LNGTHO/DFLOaT(JMAXI) 

DO  55  J»I,JMAXl 
LNGTH(J)»CFLOAT(J).DLNGTH 
55  LHLSQ(J)«(0.5D0*LNGTH(J))«P2 

N9C0LL  » .NOT.  COLL 
NOLOSS  p .NOT.  LOSS 
NOBURN  ■ .NOT.  BURN 

DATA  E,  BOLTS,  C / 4.8032  D-IO,  J.6022  D-B,  2.0970  D.10  / 
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■ 


OATA  MO,  ME,  MT  / 3.3453  0-24,  4.1045  0-28,  5,0179  0*24  / 
OAT*  ZO,  ZT  / I.OO,  2.00  / 

OAT*  gamma  / 1,666600067  DO  / 

OAT*  OPS  1 2 / 0.00  / 

DATA  TEL8SS,  NOLOSS,  NTL08S  / 0.00,  0.00,  0.00  / 

DATA  8,  RTEST  / 1.0-6,  1,0-6  / 

OAT*  At,  A2 , AJ,  AO.  45  / 2,621  01,  1.54  07,  1.46  0.6, 

1 1.247  00,  6.47  D5  / 

OATA  Z / 0.00  / 

Pi  ■ 4.00*DATAN(1.00) 

S04TPI  • DSQPTTPI) 

RT2PI  ■ OSOPT (2,00*PI ) 

TNOPI  • 2.00‘PI 
POURP I » #,00»P1 

eictp|68.doap: 

c 

C COMPILE  TABLE  OP  MISCELLANEOUS  CONSTANTS, 

c 

0M 1 ■ GAMMA  • i.OO 

KBETA  ■ 0S0RT(B8LTZ/(GMl6HE«C«C>> 

XBREM  ■ 5.350-24 

KBUPN  ■ 18.14D3*B8LTZ*0H| 

KCNVPT  » C/1.01 
MCE  ■ I ,0-7»M0 
KCVCL  ■ 3.480-16 
K0*0 .2500 
RE  « B8LTZ/GM1 
RP»1.0T/(TMOPI«LNCTMO) 

KPLLI*  ■ 1,08/PJ 

KID  » E»Z0/(TM8PI»KCNV»T) 

KIE«-E/{Tn6P1»KCNVPT) 

KIT  » E*ZT/(TNOPI»KCNVPT) 

KKIN»PI*«H8aLNGTH0/2.D0 
KMU  • P BURP] » 1 .0-4 
MU  ■ P8UPPI *1 ,0-4 
MNEUTR  » 0.00 
KNKT  ■ 1.0-7»B8LTZ 
KP  • 1.00/16. 00*PI> 

KPm I ■ 1.  0-e#TWPPl«ME*C/E 
KPP8P  « mq«c/(ZD*E) 

KP  ■ 0SQPT(3.00*HD*B8LTZ)*C/E 

KRE8»1 ,02«P0URPI 

KP6T»PI«Ph8«8MEGA**2*LNGTM0 

KT  ■ CM1/88LTZ 

KTD  » M0*GM1/B8LTZ 

KTE  ■ me*GM1/08LTZ 

KTT  « MT«Gmi/B8LTZ 

KVD  ■ 1 .08»Z0«E/(Tw0PI*MO*C) 

K VT  ■ 1 ,D8«ZT»E/(TW8PI*MT*C) 

KV8L  ■ TM8Pl«PI 
KYIELD  ■ KBUPN/G“1 


c 

c CALCULATE  CONSTANTS  TO  BE  USED  IN  COLLISION  RATES. 

e 

LOGOEO  ■ -OLOG(DSORT(FOURPI/B8LTZ,*3)»ZC«E««3J 

LOGOTO  • -DL8G(DS0RT(F8URPI/B0LTZ)«ZD*ZTaE»»3*(M04MT)/(M0«MT)> 

L9GTE0  ■ -OLOG(DSORT(E8URPI/B8LTZa»3)*ZT«E»*S) 

e 

CNUSDE  ■ (4.06*PT2PI/3.00)*ZO**2*E**4*(1.00  ♦ M0/ME)*ME*«1 .5/ 

1 (M0*«2aB8LTZ**1.5) 

CNUSOT  ■ 4.p6*PI*(ZD*ZT*E**2)«*2*(i.Do  ♦ M0/MT)/M0**2 
CNUSTE  ■ (4.Do*PT2PI/3.00)*ZT**2*E**4*(I.OO  ♦ MT/ME)*ME»«1.5/ 

1 (MT**2*B6LTZ**1.S) 

c 

CNUTOE  ■ (B, 00 *PT2P 1/3.00 )*(Z0*E*E)**2*0SGPT (ME)/(Mo,60LTZ«« 

CinkSV  1.5D0) 
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ot,u‘u;'CNl)T0T  • FfluRPI*IZD«ZT»EB£)B*2BGNl/{SCLTZ*ND) 

CNUTTC  ■ (e.D0«RT2PI/3.00)*(ZT*e*E)**2*CSQRT(He)/(MT*eet.TZ** 
i i.oooj 

e 

C OETeRMJNE  LINER  CONSTANTS. 

c 

Rioso  ■ R 1 0«R  1 0 
R20S0  ■ R2C»R20 
ROSO  ■ R20SQ  - R10SQ 
RhSQ  * RK»R* 

NTCAP»SPHT*PI*R0SQ*LNGTHO»Rh0 
PS  120  ■ 320* (RN8Q  • R20SQJ/KFLUX 

c 

c determine  initial  ring  and  plasma  parameters'. 

c 

VDO  ■ -DSGRT(2'.DO*BOLTZ*MOO/MC) 

Cl  « PI»Rmajor/1  .02 
C2  « 5,D0*PPO 

C3  * 88LTZ* (TOO  ♦ NRATIOtTTO  ♦ (ZD  ♦ NRAT 1 8 »ZT ) *TEO ) 

C«  • (RID  ♦ KIT*nRaTI8.VRaTI8)*(I.D0  - IRATISJbVDO/RMAJO 

IRING  ■ C3/(Cl*C4) 

Rmjno  » -IRING/C2 
RMIN8R  « RHINO 

V8LUME  ■ KV(JL*R«AJ9R*RMJN8P**2 
VGMJ  * V8LUME«*G"1 
NDTBTO  * I R INC /C  4 
NDO  ■ NDT8T0/V8LUME 
NTT8T0  ■ NRATI8*NDT8TC 
NTO  ■ NTTOTfl/VPLUHE 
VTO  » VRATI8»VD0 
IDO  • KID*NDT8TO*VDO/RMAJO 
J TO  ■ KIT*NTT8TO*VTO/RMAJO 
IEO  ■ IRING  - IDO  - ITO 
RJO  ■ RIO 
C 

C THE  INITIAL  VALUES  OF  THE  LINES  CURRENTS  AND  NKT  ARE  N8H  FOUND 

c through  a series  of  algebraic  equations  designed  to  insure 

c that  Zero  flux  thrfaos  perpendicular  to  the  liner  at 

C SPECIFIED  POINTS. 

c 

JMl«jMAX-2 
DO  41  J* I , JM  A X 
DO  41  JJ»1,JMAX 
«1  ZA (J, JJ) «0 . 00 

DO  44  Jb 1 « JM 1 

ZA(J,J)»KMU«OLflG((LNGTM(JBt )bLNCTH(J))/(LNGTH(Jb1) 

1 -LNGTH(J)))/(TH0PIbLNOTH(J)) 

IF(J.EQ.t)  GO  TO  44 
JPbJ-1 

DO  43  JJ*1,JP 

D ■ (LNGTH ( J* 1 1 - LNGTH(JJ)J  / 2.00 
F ■ LNGTH(JJ) 

43  ZA(J.JJ)  « (KMU/(TW0PI*LNCTH(JJ)))  * DLOG{(DbF)  / D) 

44  CONTINUE 

DO  51  J«1,JMAX1 

ZA(J»AX, J JbOPMLRS (RHINO, RH A jo, RJO, LNGTH (J).LHLSQ(J)) 

51  ZAfJHAXt , J)«DAMLRS(RMINO,RMAJO,RJO,LNGTH(J),LHLSO(J)) 
ZA(JMaX1,JHAX)b2,OO/(RMIN0*IRING) 

ZA(JMAX,JH*x).1 ,D0/(RNAJ0bIRING) 

DO  52  J»1,JM1 

OENeM«(LNGTH(jin/2.CO)*B2*(RJ0-RMAJ0)»*2 

ZR(J)«.(KMU*IRiN0/F8URPI)BLNGTHCJM)/DENeM 

52  CONTINUE 

ZR(JMAX1 )«-o,5DO*IRINGbDALR(RMINO,RMAJO) 

ZR(JH*x;«.0.5C0*IR!NGbDBLR <RMINo, Rma JO >-KCPbNDTOTObVDO***/ 

i (rhajo*iringi 
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JMAX2«JHAX»*2 

CALL  DGELGCZR,ZA,JHAX, JMAX2,1,J,D-6,JIER> 

PRINT  42,  J1ER 

it  forhatc  • # i x# 1 ierrgr*  1 » it) 

PRINT  46,  (ZR ( J) « J* 1 » JHAX ) 

«6  FORMAT ( IPSE  12,4) 

4T  CONTINUE 

DO  S3  J«t,JHAXI 
SI  ILlNERIJIiZRU) 

NKT0*ZR (JHAX) 

c 

PHIO'L* (RHINO, RHAJ0)»I»:N6 
DO  54  J* I , JHAX t 
D»LNGTH(J) 

F»LHLSO(J) 

Phi 0»PHI0«HlRJ (RHINO, RHAJ0,RJ0,O,n«lLINER(J) 

PSIl  (J)»LL(RJO.O)*ILINER(J)*HLR»{P.HINo,RHAJO,RJO,0,n*I«INft 
DO  62  JJ*1,JMAX1 
lP(J.EO.JJ)  GO  TO  62 

PSIKJ)  ■ PSH(J)  ♦ HLL(RjO,LNOTH(J),LN6TH(JJ))»ILlNERfJJ) 
62  CONTINUE 

V(J*IR)»PSII(J) 

54  CONTINUE 

NUO  ■ CNUSDT*2.DI«NT0/(0AB3(VC0  - VT0))**3 
UiO  • RHAJOR4NUO/PCA 
IF  (PORMRD)  UlO  • *UlO 
C 

C INITIALIZE  THE  DEPENDENT  VARIABLE  (VI  ARRAY. 

C 

Y(l)  ■ TIHE 

Y C 2 ) • R1080 
Y(l)  ■ R I 0 *U 1 0 
Y(4)  « 0.D0 

Y (5)  ■ PSI20 
Y(6)  ■ 0.00 

Y (7 ) ■ RH*Jo«VDO 
Y (8 ) « RHAJo«VTO 
Y (9)  ■ TC0»VGHl 
Y ( 1 0 ) ■ TEO«VGHl 
Y ( 1 1 ) » TT0«VGM1 
Y ( 12)  ■ NDTOTO 
Y (13)  • NTTOTo 
Y(14)  ■ 0.00 
Y ( 1 5 ) ■ O.Dc 

Y c 16)  ■ 0.00 
Y (I  7)  ■ O.DO 
Y ( 1 8 ) « O.Oo 

Y ( 19  ) ■ 0.00 

00  60  J*I,JHAxi 
*0  Y ( 19aJ)«P$I I ( J) 

Y(JT0TAL)»PHI0 

c 

C COPY  INITIAL  VALUES  INTO  SAVE  ARRAY. 

C 

DO  SO  Jat , JTOTAL 
50  YO (J)  • Y (J) 

C 

C INITIALIZE  EXTRAPOLATION  PROCEDURE  FBR  FIRST  0UESSE8  USED  IN 
C NE"T0N«RAPHS0N  ROUTINE  ON  EVEN  STEPS  OF  R»K>C. 

C 

DO  61  J ■ 1,  4 
RHAJl(J)  ■ RMAJO 

Rha J2 ( J ) ■ rhajo 

RHINKJ)  ■ RHINO 
*1  RHIN2CJ)  ■ RHINO 

t ) n<?r:v  otne*  ■ dt 
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- - ftETUR(g 
e 

ENTRY  08PLAYCSNAME) 

DATA  |V,  SYO,  10V,  *11,  (RING,  SP|,A8,  SNU,  SQNU,  I0ATA  / 

1 <Y  <,  • YO  *,  'OV  »,  'Ll  »,  'RING*,  'PLAS*.  'NU  », 

2 'CNU  <,  'DATA*  / 

OATA  SBLANK  / • ' / 

00  100  J ■ 1,  33 

100  HINE(J)  ■ SBLANK 

c 

IP  (SNAME  .EG.  SY  1 GO  TO  110 

IP  (SNAME  .EG.  SYO  ) GO  TO  120 

IP  (SNAME  .EG.  SDY  ) GO  TO  130 

IP  (SNAME  .EG.  SLI  ) GO  TO  140 

IP  (SNAME  .EG,  SRING)  GO  TO  ISO 
IP  (SNAME  .EG.  SPLAS)  GO  TO  1*0 
IP  (SNAME  .EG.  SNU  ) GO  TO  170 

IP  (SNAME  .EQ.  SCNU  ) GO  TO  100 

IP  (SNAME  .EQ.  SDATA)  GO  TO  190 
PRINT  101,  SNAME 

101  FORMAT ( *0 1 « 'OSPLAY  CALLED  WITH  INVALID  ARGUMENT  ',  AO  //) 
RETURN 

C 

C DEPENDENT  VARIABLES. 

C 

110  PRINT  111 

111  PORMATd-  Yd)  Y (2)  Y (3)  V (5 ) Y(6)  Y(7)  ', 

1 ' Y (9  ) Y (10)  Y ( 11 ) Y (12)  Y ( 13)  YtlO)  Y(l5)', 

2 ' Y ( 1 A ) Y ( 1 7 ) V ( 16 ) ' ) 

CALL  SCNVR5(SLINE,_  Yd),  Y(2),  Y (3) , Y(5),  Y(»),  Y (73 , Y(9), 

1 Y (1 0)  , Y(ll),  Y ( 12) , Y(13) , Y ( 1 0 ) , Y(l5),  Y(16>,  Y(l7), 

2 7(1833 
PRINT  199,  SLINE 
RETURN 

C . . 

C INITIAL  VALUES  OP  DEPENDENT  VARIABLES. 

e 

120  PRINT  121 

121  PORMATd-  YO(t)  YO  (2)  YO  (3)  Y0(5)  Y0(6)  Y0(7)  ', 

1 'VO  ( 9 j YO  (10)  YO(ll)  YO  (12)  YO  (l3>  YO(lfl)  Y0(l3)', 

2 ' YO (1 6}  YO ( 1 73  YO  (1 8 3 ' } 

CALL  8CNVR5 (SLJNE , YO(l),  Y0(2),  YO(3)4  YO  (53 , YO(A),  YO (73  # 

1 YO  (9 ) , YO(tO),  YOdl),  YO  ( 1 2) , Y0(l3),  YO(lO),  Y0(l5), 

2 Y o ( 1 6 ) , YO (17) , YO (IB) ) 

PRINT  199,  SLINE 

RETURN 

C 

C DERIVATIVES  op  DEPENDENT  VARIABLES. 

c 

ISO  PRINT  131 

111  PORMATd-  DY  (I ) DY  (2}  0Y(3)  DY  (5  3 0Y(*)  OY  (7)  ', 

j ' DY  (9  3 DY  (1 0 ) OY(U)  DY(12j  DY(l3)  DY(Ifl)  OY  ( 1 53 ' , 

2 « DY (1  A)  DY  (17)  DY (18)') 

CALL  SCNVRS (SLINE,  0Y(1),  0Y(2),  0Y(3),  DY(5),  DY(6),  OY (7 3 , 

1 DY (9 ) , OY (| 03 , DV ( 11 ) , DY(12),  0Y(13),  DY(IS),  OY(lS), 

2 DY (1 6) , DY (1 7) , DY (18) 3 
PRINT  199,  SLINE 

RETURN 

e 

C LINER  PARAMETERS  AND  ENERGETICS. 

c 

ISO  PRINT  1 « 1 

HI  PORMATd-  R 1 R2  U|  U2  G1  G2  ', 

1 • B1  82  ETA  P8I2  RSKIN  PSIl(l)  PSIK2)', 

2 ' PStl  (3)  PSIl(S)  PSIl(Sl') 
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1 


1 *812,  R J , PSIl(l),  Pill  (2) , *81 1(3),  *811  ( * ) * *811(5)) 

PRINT  HI,  ILINE 
PRINT  142 

142  PORMAT ( 1 " CTOT  £«IN  ER8T  E*0T  £HAG  EMAG2  '• 

I 'ERIN6  EBEAH  EDDIR  EOTH  E*LA8  EETH  ettm 

1 ' ETOIR  ENET  O') 

CALL  SCNVRSdLINE,  ETOT,  ERIN,  £R8T,  E*8T,  £MAC,  EMAG2  > 

1 CRING,  E*EAR,  C00I*#  EOTH,  E*L*8,  EETH,  ETTM,  ETOI*# 

1 CNET , 0) 

PRINT  199,  ILINE 
PRINT  143 

1«1  P8RNATC-  EH  AS  ( 1 ) CHACC2)  EMAG(3)  CHAG(4)  CMAG(S)  XL  1 •# 

1 ' X L 2 1L3  IL«  IL5  RLt  *L2  «LS  '# 

2 ' RL«  RL5  TU • ) 

CALL  SCNVRSdLINE,  CMAG1  ( I ) « £HAG l (2) « CHAG1 (3) # CHAGl(«)> 

1 EHAG1 (5 ) , ILINER(l)#  ILIN£R(2),  2,2,2, 

2 RLINER(t).  RLINERC2),  RLIHER (3) # RLINERd), 

I «LINER(5>,  TL) 

PRINT  199,  ILINE 
RETURN 
C 

C PARAHETERS  8*  I8N  AND  IMAGE  RINGS. 

c 

150  PRINT  151 

151  P8RMAT ( 1 • RmajoR  RNIN8R  IRING  10  L*ING  MLR 1 '# 

1 ' PHI  PD  VD  HD  NDT8T  ND  TO  '# 

1 'V8LUMC  ITER  RATE') 

CALL  SCNVRSdLINE,  RMAJOR#  RMINBR,  IRING,  10#  LRING,  MIR1I, 

1 PHI,  PD,  VD,  HO,  NDTOT , NO,  TD,  VOLUME , ITER,  RATE) 

PRINT  199,  ILINE 
RETURN 
C 

c target  plasma  and  radiation  parameters. 
c 

1*0  PRINT  141 


1*1 

P8RMATC-  NETpT 

NT  TOT  VOLUME 

NE 

NT 

HT  ', 

1 

•VE  VT 

PT0T  TE 

TT 

PE 

PT  ', 

2 

• IE  IT 

TIELD') 

CALL  SCNVRSdLINE, 

NET8T,  NTTOT, 

VOLUME, 

NE,  NT, 

HT,  VE,  VT, 

1 

PT8T,  TE,  TT, 

PE,  PT,  IE,  IT, 

UELO) 

PRINT  199,  ILINE 

PRINT  162 

1*2 

PORMaT ( • • beta 

BP  PCTCL 

PBREM 

PRAO 

NURAD  ', 

1 ' ERAO  OTERAD  PP8L  RNUBTU  NC8UNT  BETaPL  E8MMIC  ', 

2 ' KAPPA  PVGPRD  8URNUP ' ) 

CALL  SCNVRSdLINE,  BETA,  B*,  *CTCL,  *6*1*,  PRAO,  NURAD , ERAO, 
t OTERAO,  PP8L,  RNUBTU,  NC8UNT,  BETAPL,  E 6HMIC , KAPPA, 

2 PVGPRD,  BURNUP) 

PRINT  199,  ILINE 
RETURN 
C 

C TRANSPORT  RATE!. 

C 

170  PRINT  171 

171  PSRMATC-  NUSOE  nusot  nuste  NUSED  NUITD  NU8ET  ', 

1 ' VOT  VOE  VET  CLOG  DHCPE  OHD/T  DMT/e  ', 

2 ' DHE/O  DMT/D  OHE/T') 

CALL  SCNVRSdLINE,  NU9DE , NUSOT,  NUSTE,  NUSEO,  NUSTD,  NuSET, 

1 VOT,  VOE,  VET,  EL8G,  DHDSE,  DHDST,  dhTIe,  DHElD,  DMTSO, 

2 OnEIT) 

PRINT  199,  ILINE 
PRINT  172 

172  P8RMAT ( ' • NUTOE  NUTDT  NUTCT  NUTED  NUTTD  NUTTC  ', 

1 'OTO/E  OTO/T  OTE/T  OTE/D  DH/D  DTT/E  LBCLDE  ' , 

2 ' LBGLOT  LOGLTE') 

■rori'  fl-T  CALL  SCNVRSdLINE,  NUTOE,  NuTOT,  NUTET,  NyTffr,  NUTTD,  NuTTE, 
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C 

C 

c 

1 7J 


175 


C 

e 

c 


150 

m 


i*i 


i*i 


t** 


c 

l«« 


1 “ r DTD  JE , 0Td*T»  dtejt,  dtesd.  dttjc,  cttie,  logloe.  loglot* 

2 LOGLTE.  Z) 

PRINT  19*,  JUNE 

INDUCTANCES. 

PRINT  l 73 


FORMAT!'-  LLl 

LL2 

LL3 

LL«  LL5  THIS  SPACE’ 

1 

' FOR  RENT  ') 

CALL  SCNVR5!JlInE, 

LLl, 

L L 2 a l.  L 3 $ 

LL®,  Ll3,  Z,  Z , Z,  Z,  It 

i 

Z,  Z,  Z,  z.  z, 

ZJ 

PRINT  19*,  JLINE 
PRINT  175 

FORMAT!'.  LL5 

ML  12 

MU3 

ML  J R ML  1 5 ML23  ' 

F0HM»TC'»  LL5  "'-It  “I.1J  MLJ«  ML15  ML2  3 

1 'ML2«  ml 25  ML3«  ML3S  MLR5  MLfiJl  HLRS2  '* 


2 ' MLRJ3  MLR  JR  M|_R  J5  ') 

CALL  SCNVRStfLlME,  LI.  3 , ML12,  “1.13,  “LI*.  ML15,  ML2J.  “l2«, 

1 ML25,  ml 3R « Mi.35,  mlr 5 , MLRJ’.,  MLRJ2.  mlRJJ,  MLRJtt,  MLRJ5) 

PRINT  199,  JLINE 

return 

constants  used  in  evaluating  transport  and  radiation  rates. 

PRINT  181 

FORMAT!'-  L8GDEO  LOGCTO  LCCT20  KTD  KTE  KTT  ', 

1 'CNUSOE  CNUSDT  CNUSTE  CNUTDF  CNUTOT  cnutte  rbeta 

2 ' K8REM  KCTCL  NT’) 

CALL  SCNVR5C JLINE,  L OGCEO , LSGDTO,  LOGTEC*  KTD,  KTE,  KTT, 

1 CNUSCE,  CNUSDT,  CNUSTE,  CNUTDE , CNUTCT , CNUTTE,  KBETA, 

2 K0REM,  KCTCL , KT> 

PRINT  199,  JLINE 
RETURN 

PHYSICAL  CONSTANTS  ANC  ST-SEC  COMBINATIONS  THEREOF. 

PRINT  191 

FORMAT!'-  £ S9LTZ  C MD  me  MT  ', 

1 ' ZO  ZT  P4  LENGTH  TELOSS  NDLOSS  NTLOSS', 

2 ' CPS  1 2 OMEGA  GAMMA'} 

CALL  5CNVR3( JlINE.  E,  30LTZ,  C,  "D,  ME,  MT,  ZD,  ZT,  RN, 

1 LNGTho,  TELOSS,  NDLOSS,  NTLOSS,  CPS  12 , OMEGA,  GAMMA) 

PRINT  J99,  JUNE 
PRINT  192 

FORMAT!'-  KBUPN  KCNVRT  KCF  KD  KE  Kr  ', 

1 'KFLUX  me  K IE  KIT  KKJN  KHU  KNEUTR  ', 

2 ' KNKT  KP  KPPOP') 

CALL  SCNVR5(»LINE,  K8URN,  KcNV»T,  KCF,  KD,  KE,  KF,  KFLUX, 

1 Kic,  kIE,  KIT,  KK I N , KHU,  KNEUTR,  KNkT,  KP,  KPROP) 

PRINT  |99,  JUNE 
PRINT  193 

FORMAT!'-  KPHI  KR  K»ES  KROT  KVD  KVT  ', 

1 ' *VC..  KVD  KVT  KVIELD  CR  RTEST  SPHT  ', 

I ' HTCAF') 

CALL  SCnvR5(Jl1UE,  KpHI,  K»,  K»ES,  KROT,  KVp,  KVT,  KVOL,  KVD, 

I KvT,  WIELD,  D,  "TEST,  SPHT,  MTC A P , Z,  Z) 

PRINT  '.99,  JLINE 
PRINT  1911 

PeR“ATP-  DT  D’PLOT  DTFILM  DTCLMP  TLAST  T 

1 'TPL8T  TFILM  TDUMP  IRATie  VRATIO  NRATIO') 

call  BCNVRSCIlINE,  DT,  DTPLOT,  DTFILM,  CTDUMP,  TLAST,  TIME, 

1 TPlOT,  TFILM,  TDUMP,  I RA T 1 0 , VRATIO,  NRATIO,  Z,  Z,  Z,  Z> 

PRINT  199,  JLINE 
RETURN 
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ENTRY  8UTPUT 
RETURN 
END 

C 

RE*L  FUNCTION  FA*8(A,B,C,NKT) 

IMPLICIT  RfAL*8(A-I,K-Z),  INTE6ER*«CJ) 

DIMENSION  YAH, 3) 

COMMON  / INDEX  / JMAX,  JMAXJ 
COMMON  / PBLOCK  / P 
COMMON/  RADII/  RMJN8R,RMAJ0R,RJ 
COMMON/  LENGTH/  LNGTHC5),  LHL80<5) 

COMMON/  ARRAY/  Y <25  3 * OYC2S),  Y0I25) 

COMMON  / INCO  / MU,  KMU 
common/  F lux/  PHI,PSIl(S) 

COMMON  / CURREN  / R(3) 

LR ( A , 8 ) » KMU*8*(DLOG(O.DO«8/A)  - 2.00  ♦ KD) 

DALR(A.B)  ■ -kmu»b/a 

DBLR(A,B)  » mmu«(DL0G(8.D0»B/A)  - 1.C0  ♦ KD) 
LL(C,D)«MU«C»(DL0G(16.DO»C/D)  - 0.5DC) 
DCLL(C,D)»MU*tDLOGtib.OO*C/D)  ♦ 0.500) 
LN(D,F)*0.5D0»cD/F)oDLO&CD) 

FAC <D,F)«{((F-D)**2)/(2.D0*DaF))«DL0G (CABS (F-D)  ) 
MLL(C,D,F)PMU*C»(0L8G(16.00«C)-LN(D,F)-LNCF,D)4FAC  fD.F). 0.500) 
OCMLLCC.O, F)»MU«COL&GC16.DO»C)-LN(r),F).lN(F,0)*FAC(P,Flf0.5DO) 
MLR(B,C,F)»MU«D8QRT (B»C)«(.5D0»DL0G(6.«C1 *6.C/( (B-C)«»2  ♦ 

1 Fj)  - 1 .00  - {(B-C)/OSORT(F))*rATAN(CSCRTCF)/(B-C))) 

OBMLR(B,C,F)i,5DO»MU»DSQRT(C/B)*(.5DO»DLOG(6.R01«B»C/t(B*C)»*2 
l ♦FJJTC(C”J,.DO*B)/OSQRT(FJJ*DATAN(05CRT(F)/(B*CJJ) 

0C^LRtB,C,E)«.5D0»MLI«0SQRT(B/C)»(.5D0*DLeG(b.ttDt *B»C/((B-CJ»p2 
I ♦F))*((J,00*C*B)/OSQRT(F))*DATAN(DSCRT(F)/(B"C))) 

ML«*(A,B,C,0,F)«MLRCB,C,F)A(DSORT(LR(A,e)*LLtC,D))-MLR(B,C,F)) 

1 /(1.D0*((C»B)/A)**P) 

DANLRSCA,8,C,0,F)«(.5D0*DALR <A,B)*DSORT<LL tC,D)/LR(A,BJ)T(p/A) 

1 *C(C"B)/Aj»*P*(DSQRT(LRtA,B)«LLfC,0)J»MLR(B,C,F))/fl.DOA 

2 C(C-B)/A)»*P))/(i.DO  ♦ ((C-B;/A)«.P) 
0BMLRJ(A,B,C,D.F)«OBMLR<8,C,F)a{.5O0»OBLR(A,B)*08QRT(LL(C,D)/ 

1 LfifA.B))  ♦ CP/(C-B))*((C-B)/A).«P»fCSQRT(LRCA,B3«LLtC,0})- 

2 MLR  (6 , C ,F  } ) / ( 1 '.DO  ♦ (<C-B)/A)**Pl)/(l.DO  ♦ ( <C«B 3 /A ) . «P > 
OCMLR8(A,B.C,O.F)rDCMLR(B,C,Fj»(.5DO*OCLL{C,D)*OSQRT{LR(A,B)/ 

1 LLCC,0))-(P/j[C»B))»CCC-B)/A).«P»(D8CRT(LR(A,B  j»LL(C,D))-  • 

2 MLR(B,C,F))/{1.D0  ♦ ({C-B)/A)*.pj)/(I,00  ♦ ( (C-BJ/A ) **f) 

C A SERIES  OF  ALGEBRAIC  EQUATIONS  IS  SOLVED  TO  FIND  THE  CURRENTS 
C IN  THE  LINER  AND  RING. 

KD  ■ .2800 
MINR«A 
MAJR«8 
LIR*C 
LIR2«LIR 
JMAXMiJMAX-l 
DO  1 J»1,JMAXM 
l R ( J ) ■ V C Up  J) 

R (JMAX)rY (IRpJMAX) 

YA(1,1)«  LL(LIR,LNGTM(t)J 
DO  2 Jt»2,jMAXM 
A1*LNGTH(J1) 

YA(Jl,Jn»LL(LIR#AM 
JIMI»JI-l 
DO  2 J2«t,JIMl 
F»LNGTH(J2) 

YAfJl,J2)»MLL(LIR,A1,F) 

YA(J2,Jl)«YA(JI,J2) 

I CONTINUE 

DO  « J»1,JMAXM 
AlPLNGTH(J) 

OpLHLSOCJ) 
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i qn l?3Ty*um*)‘* j5»mlr»<min«,majr,l:r,ai,d) 

<T  YA(J, JMAX)*YA(JMAX, J) 

YACJHAX,  JMAX)bi.R{HInR,NAJR) 

JMAX2bJMAX««2 

CALC  OGELG  (R,  YA,  JMAX,  JM*X2,  1,  1.0-6,  JIER) 

FA  b6.5D0*R<JNAX)«*2*DALRIA,B)*2.D0*NKT/A 

09  20  Jat.jMAXl 

OaLNGTH(J) 

FbLHLSQCJ) 

20  FA  aFA  ♦R(JMAX).R{J)*0AMLR»(A,B,C,0,F) 

return 

END 

REAL  FUNCTION  FB*8(A,B,C,NKT) 

IMPLICIT  REAL*«(A-I,K-Z),  INTEGER*fl(J) 

DIMENSION  YACJ.J) 

COMMON  / INCO  / HU,  KMU 
COMMON/  COEFF/  KCF.NOTOT,  RVO 
COMMON  / PBLOCK  / R 
COMMON  / CURREN  / RC3) 

COMMON/  LENGTH/  LNGTH(5),  LHLSQ(S) 

COMMON  / INDEX  / JMAX,  JMAX 1 
COMMON/  array/  Y (25 J t DY  f 25 ) , Y0C25) 

COMMON/  RADII/  RMJNOR,RMAJOR,RJ 
COMMON/  FLUX/  PHI.PSIl C5) 

LR  ( A , 8 ) a KmU*0*(DLOG(8,DO*B/A)  - 2. CO  •*  KD ) 

DALR  ( A , B ) « -KMU-B/A 

DBL  R ( A , B ) a KMU* ( DL OG ( 8 . DO *0 /A ) - l.CO  ♦ KD) 
LL(C,D)*MU*C*(DLOG(16.D0*C/D)  • 0.500) 
DCLLCC,D)*mU*(OLOGU6.DO»C/D)  ♦ 0.500  ) 
LNCO»F)»0,5DO*(D/F)«DLOG(D) 

FAC(D,F)*({(F-D)**2)/(2.DO*0«F))«DLOG(OA8S(F-0)) 
MLLtC,D,F)aMU.C*fDL0G(16.D0»C)-LN(D,F)-LNCF,O)*FAC(O.F)-O.5DO) 
DCMLL(C,D,F)a-U»{DLOG(l6.DO*C)-LN(D,F).LNCF,0)*FAC(D,F)*0.5DO) 
MLR(B,C,F)aMU«DSORT(B«C)*(.500»DLOG(6,OCl*8*C/((B*C)»*2  ♦ 

1 F))  - 1.00  - ((B-C)/OS0RT(F)i*pATAN(CS0RT(F)/(B-C))) 

DBMLR(B,C,F)a>jD0»MU*DS0RT(C/B)*f.5DC*0L0G(6.RDl«B«C/((B-C)»*2 
l ♦F))A((C-3.OO*B)/OSORT(F)).DATAN(DSCPT(F)/(0-C))) 

DC’4LR(B,C,Fj»t5D0*MU*DS0RT(B/C)*(.5D0*DL0G(6,R01*B«C/(<B-C)*«2 
1 ♦F))f(f3.00*C-B)/DSQRT(F))*DATANtDSCRTtF)/(B-C))) 

MLR|{A,B,C,D,F)«MLR{B,C,F)*(DS0RT(LR(A,e)«LL(C,D))-MLR(8,C.F)) 
1 /CI.DO*CCC-B)/A).*P) 

DAMLRj(A,B,C,D,F)»(.5D0aDALR(A,B)*DSCRT(LL(C,D)/LR(A,8))*(R/A) 

1 *((C-8)/A)**P.cnS0RT(LR(A,B)*LLfC,0))«MLR(B,C,F))/(I,00* 

2 C(C*8)/A)**P))/(1.D0  ♦ ( (C*B)/A ) *»R) 
O8MLR|(A,8,C.D,F)aoBMLR(B,C,F)*(.5D0»DBLR(A,B)»DS0RT(LLCC,D)/ 

1 LR ( A , B ) ) ♦ tpK(C-B))»{(C-B)/A)««P.{CSGRT(LRtA,B)*LL(C,0))- 

2 Ml.R(B,C,F))/(t.DO  ♦ ((C-B)/A)*«Pj)/(l,00  ♦ ( (C-B  ) /A  ) *.P  ) 
DCMLRj(A,8,C,D,F)Rf)CHtR(B,C,F)A(.  500*DCLL(C,D)*D80RT(LR(A,B)/ 

1 LLfC,D))-(R/(C-8))*({C-B)/A).«R.(08CRT(LR(A,B  )«LL(C,D))- 

2 MLR(B,C,F))/fl'.DO  ♦ ((C*B)/A)**P))/(l.OO  ♦ ( £C-8  ) /A  ) *«? ) 

C A SERIES  OF  ALGEBRAIC  EQUATIONS  18  SOLVED  TO  FIND  THE  CURRENTS 
C IN  THE  LINER  ANO  RINO. 

KO  a , 25D0 
MINRaA 
HAjRafl 
LIR»C 
LIR2«LIR 
JMAXMbJMAX-1 
DO  1 Jal.JMAXM 
i R(j)«vaR*j) 

P(JMAX)aY(lR*JMAX) 

V A C 1 . 1 ) ■ LL(LIR,LNGTH(D) 

DO  2 J1*2,JMAXM 
AlaLNCTMIJl ) 

TA(Jl,JI)»LL(LtR*Al) 

J I “ I ■ J I • I 
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PiLNRTHC J2) 

YAtJl.J2)»MLL(LIR.M#F) 

YA(J2,J1>»YA(J1,J2) 

2 CO^TI NUE 

08  fl  J«1,J*AXN 
*l»LNGTH(J) 

0»LHLS0(J) 

Y*(JM*X,J)»MLBJ(MINR,MAJR,l:R, Al,D) 

a ya(J,jmax)»ya(jmax,J> 

YAfJMAX, JM*X)«LR(M1NR,MAJR) 

JMAX2> JHAX**2 

CALL  OGELC  (R,  YA,  JHAX,  JMAX2 1 1,  1.0-6,  JIER) 
FB«0.500»RCJMAX)*»2*08LR<A,B)*NKT/BtKCF*NDTOT*RVD»«2/B*»J 
00  20  J«l,JMAXl 
0»LN8TH(J) 

F»LHL$0(J) 

20  FB  »FB  ♦R(JMAX)*RCJ)oOBMLR»<A,B,C,D,F) 

RETURN 
END 
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EQUIVALENT  CIRCUIT 


Fig.  2 - Electrical  circuit  equivalent  to  Fig.  1 . A homopolar  generator 
is  used  to  energize  the  driving  coil. 


37 


BOOK,  TURCHI,  AND  STEIN 


38 


NRL  MEMORANDUM  REPORT  3827 


39 


NRL  MEMORANDUM  REPORT  3827 


41 


ENERGY  CMJ)=o 


BOOK,  TURCHI,  AND  STEIN 


I ig.  7 Magnetic- energy  [from  Eq.  (9 ) J , liner  kinetic  energy,  ion 
streaming  energy,  and  total  particle  thermal  energy  vs.  t. 
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